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I.    Introduction 

This  report  contains  user  instructions,  a  listing  and 
documentation  for  a  microcomputer  BASIC  program  that  can  be  used 
to  compute  an  estimate  of  the  probability  that  a  magnetometer 
based  detection  system  such  as  the  AN/ASQ-81  will  detect  a 
submarine  magnetic  dipole  field  during  an  encounter. 

The  program  generates  detection  probabilities  based  on  two 
encounter  models.   In  the  first  encounter  model,  the  detection 
system  uses  a  crosscorrelation  detector.   In  the  second 
encounter  model,  the  detection  system  uses  a  square  law  detector, 
Relative  to  operationally  realizable  values,  probabilities  based 
on  the  first  model  represent  upper  bounds  and  those  based  on  the 
second  model  represent  lower  bounds.   For  both  encounter  models, 
the  signal  is  proportional  to  the  square  of  the  magnitude  of  the 
projection  of  the  dipole  field  on  the  earth  magnetic  field  and 
the  magnetic  noise  does  not  change  with  changes  in  the  position 
of  the  magnetometer.   Also,  for  both  encounter  models,  an 
encounter  is  a  straight  line  encounter  with  constant  vertical 
separation. 

The  encounter  models  can  be  interpreted  as  models  of  a 
magnetic  anomaly  detection  (MAD)  system  on  an  aircraft  that  is 
moving  with  constant  course,  speed  and  altitude  in  an  encounter 
with  a  submarine  moving  with  constant  course,  speed  and  depth. 
Or,  they  can  be  interpreted  as  models  of  a  stationary  magnetic 
anomaly  detection  system  in  an  encounter  with  a  submarine  moving 
with  constant  course,  speed  and  depth. 


The  program  parameters  include:   encounter  latitude  and 
longitude,  submarine  induced  magnetic  moments,  submarine 
permanent  magnetic  moments,  submarine  course,  speed  and  altitude, 
magnetometer  course,  speed  and  depth,  encounter  lateral  range 
(the  horizontal  range  at  the  closest  point  of  approach  in  a 
straight  line  encounter)  and  false  alarm  rate.   (In  the  program, 
a  false  alarm  is  the  event  that  the  detection  system  classifies 
noise  as  a  dipole  signal.) 


II.   Program  User  Instructions 

As  listed  in  Appendix  10,  the  program  can  be  run  under  a 
BASIC  language  that  is  compatible  with  IBM  PC  BASIC.   If  the 
listing  is  used  to  enter  the  program  through  the  keyboard,  then 
the  program  should  be  saved  with  the  name  MAD.BAS  or  the  value  of 
N$  on  line  40  should  be  changed  to  the  file  name  under  which  it 
is  saved. 

The  program  contains  user  instructions  in  the  form  of  query 
and  parameter  limitation  messages.   As  an  example  of  the  former, 
after  starting  the  program  under  BASIC,  the  following  message 
should  appear: 

Magnetic  Anomaly  Detection  (MAD)  Lateral  Range  Function 

generate  or  print  a  program  data  file  (g/p)? 

By  entering   p,   data  can  be  printed  from  a  program  data  file 
that  was  generated  by  the  program.   By  entering   g,  a  program 
data  file  can  be  generated  for  a  set  of  user  specified 
conditions.   With  either  response,  a  sequence  of  additional 
queries  is  displayed.   These  queries  require  either  an  indicated 
response  or  a  parameter  value  as  the  input.   If  the  initial 
response  is   g,  the  sequence  includes  queries  whose  responses 
determine  whether  or  not  an  auxiliary  data  file  will  be  generated 
that  can  be  used  for  future  input  of  magnetic,  processing  or 
kinematic  data.   In  particular,  the  first  query  in  the  sequence 
gives  the  option  of  using  a  combined  magnetic,  processing  and 
kinematic  data  file.   If  the  response  indicates  that  it  should 
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be  used  and  the  file  is  available,  the  parameter  values  that 
remain  to  be  entered  in  order  to  generate  a  program  data  file  are 
the  following:  the  false  alarm  rate,  the  magnetic  noise,  the 
maximum  encounter  lateral  range  and  the  lateral  range  step.   The 
combined  file  should  be  used  only  if  the  effect  of  varying  just 
one  or  more  of  these  parameters  is  desired.   If  the  response 
indicates  that  the  file  will  not  be  used,  queries  concerning 
magnetic,  kinematic  and  processing  parameter  values  are 
displayed. 

After  all  of  the  program  parameter  values  have  been 
entered,  a  query  is  displayed  giving  the  option  of  generating  the 
combined  file.   Then  a  query  giving  the  option  of  generating  a 
program  data  file,  a  query  giving  the  option  of  printing  the 
encounter  parameter  values  and  a  query  giving  the  option  of 
printing  lateral  range  function  values  are  displayed.   The 
lateral  range  function  values  are  the  encounter  detection 
probabilities  indexed  by  lateral  range.   The  parameters  maximum 
lateral  range  and  lateral  range  step  determine  the  index  lateral 
ranges  of  the  encounters  for  which  probabilities  of  detection  are 
computed. 

The  program  generates  magnetic  signal  values  that  correspond 
to  points  in  time  during  an  encounter.   Following  the  lateral 
range  function  query,  a  query  is  displayed  that  gives  the  option 
of  printing  the  magnetic  signal  values  for  an  encounter.   If  the 
option  is  exercised,  the  option  is  repeated.   When  the  option  is 
not  exercised,  a  query  is  displayed  that  gives  the  option  of 


generating  or  printing  a  new  program  data  file.   If  this  option 
is  not  exercised,  the  program  ends. 

Some  suggested  guides  for  determining  parameter  values  can 
be  found  in  Section  IV  of  this  report. 


III.  Encounter  Model  Limitations 

In  the  two  encounter  models,  the  signal  is  the  unfiltered 
signal  that  is  generated  by  a  magnetic  dipole  that  moves  relative 
to  a  magnetometer.   Depending  on  the  input  filtering,  describing 
a  ship  magnetic  anomaly  as  a  dipole  anomaly  should  not  be  a 
significant  limitation  for  an  encounter  slant  range  at  the 
closest  point  of  approach  (CPA)  that  is  greater  than  one  hull 
length.   The  detection  decision  is  based  on  samples  from  a  single 
time  interval  (window)  that  is  centered  on  the  CPA.   The  length 
of  the  sampling  interval  and  the  sampling  rate  are  parameter 
values  that  are  inputs  to  the  program.   In  terms  of  signal-to- 
noise  ratio,  there  is  an  optimum  sampling  interval  length 
(integration  time)  and  sample  rate.   Although  dipole  signal 
energy  is  not  symmetrically  distributed  about  the  CPA  time,  for  a 
given  sampling  interval,  the  difference  between  the  signal  energy 
for  the  optimum  interval  location  and  the  CPA  centered  location 
should  not  be  significant  in  most  cases. 

The  encounter  model  magnetometer  noise  samples  are 
determined  by  a  gaussian  random  process.   They  are  values  of 
independent  identically  distributed  normal  random  variables.   The 
standard  deviation  of  these  random  variables  is  referred  to  as 
the  magnetic  noise  and  the  variance  is  the  noise  in  the  sense  of 
the  signal-to-noise  ratio.   The  degree  of  correspondence  between 
this  process  and  operational  noise  depends  on  the  nature  of  the 
dominant  operational  noise  sources  and  on  the  magnetometer  input 
filtering  (noise  whitening) . 


In  the  encounter  models,  the  intervals  are  adjacent  but  not 
overlapping  and  a  detection  statistic  corresponds  to  each  sample 
interval  and  its  value  is  determined  by  the  sample  values.   If 
the  value  of  the  detection  statistic  for  an  interval  equals  or 
exceeds  a  threshold  value,  a  detection  is  indicated.   The 
threshold  value  is  determined  by  the  false  alarm  probability 
which  in  turn  is  determined  by  the  false  alarm  rate  and  the 
sampling  interval  length. 

A  detector  that  used  a  moving  sample  interval  that  was 
generated  by  replacing  the  oldest  sample  by  the  newest  one  would 
correspond  more  closely  to  the  detector  in  an  operational 
detection  system.   In  a  model  of  the  detector,  since  the  sample 
windows  overlap,  the  detection  statistics  would  represent  a 
sequence  of  random  variables  that  were  correlated  over  an 
interval  equal  to  the  width  of  the  sample  interval.   Because  of 
this  dependence,  it  seems  unlikely  that  the  results  that  would  be 
obtained  with  an  encounter  model  based  on  an  overlapping  interval 
detector  would  differ  significantly  from  those  obtained  with  an 
encounter  model  based  on  a  nonoverlapping  interval  detector. 

For  encounter  lengths  of  the  order  of  a  few  nautical  miles 
or  less,  the  straight  line  encounter  condition  should  not  be  a 
significant  limitation.   In  particular,  this  should  be  the  case 
for  a  fixed  magnetometer  since,  for  a  submarine  (or  surface  ship) 
target,  vertical  separation  and  course  changes  should  be  less 
likely  to  occur. 


Other  models  are  available  that  can  be  used  as  the  basis  for 
computing  an  estimate  of  the  probability  that  a  magnetic  anomaly 
detection  system  will  detect  a  submarine  during  an  encounter. 
For  example,  one  that  is  described  in  Appendix  8  can  be  used  to 
determine  the  slant  ranges  of  straight  line  encounters  for  which 
the  detection  probability  is  equal  to  a  specified  value.   The 
parameter  values  that  are  required  to  do  this  are  an  average 
submarine  dipole  moment,  a  detection  system  capability  factor  and 
a  noise  factor.   Values  for  these  parameters  can  be  determined 
from  operational  data.   However,  the  values  are  specific  to 
averages  over  a  particular  set  of  encounter  conditions.   An 
advantage  of  the  two  encounter  models  relative  to  this  model  is 
their  adaptability  to  different  magnetic,  processing  and 
kinematic  conditions. 
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IV.   Parameter  Values 

The  magnetic  parameter  value  queries  are  generally 
explanatory  with  regard  to  the  value  that  should  be  entered. 
This  is  also  true  of  the  kinematic  parameter  value  queries. 
However,  there  is  some  ambiguity  with  respect  the  processing 
parameter  value  queries  and  the  noise  parameter  value  query.   To 
reduce  this  ambiguity,  a  brief  discussion  of  the  common 
characteristics  of  the  two  encounter  model  processing  parameters 
and  the  noise  parameter  is  given  below.   This  is  followed  by 
some  guidelines  for  choosing  these  parameter  values. 

In  both  encounter  models,  a  decision  is  made  at  the  end  of 
each  sampling  interval.   The  decision  is  either  noise  energy  was 
present  during  the  interval  or  noise  energy  plus  signal  energy 
was  present  during  the  interval.   The  sample  intervals  are 
adjacent,  equal  width,  nonoverlapping  time  intervals.   The  number 
of  samples  that  are  input  in  a  sample  interval  is  determined  by 
the  sampling  rate  and  the  interval  length. 

The  program  default  choice  for  the  sampling  rate  is  2 -MAXF 
where  MAXF  is  a  parameter  that  is  labeled  the  maximum  magnetic 
signal  frequency.   This  sample  rate  is  the  Nyquist  rate  for  an 
ideal  low  pass  filter.   However,  the  signal  in  a  sample  interval 
that  is  computed  by  the  program  represents  an  unfiltered  dipole 
signal.   This  is  a  reasonable  approximation  if  the  signal  energy 
that  is  associated  with  signal  components  greater  than  MAXF  is 
relatively  small.   As  discussed  below,  the  noise  energy  in  a 
sample  interval  should  be  considered  to  be  proportional  to  MAXF 


in  order  to  be  consistent  with  the  encounter  models.   Ideally,  a 
default  choice  for  MAXF  would  make  the  ratio  of  the  signal  energy 
to  the  noise  energy  a  maximum  for  the  sampling  interval  of  an 
encounter.   The  program  default  choice  for  MAXF  is  2 • MAXVM/MINRO 
where  MAXVM  is  a  user  estimated  maximum  encounter  relative  speed 
converted  from  knots  to  meters  per  second  and  MINRO  is  a  user 
estimated  maximum  slant  range  at  CPA  in  meters  in  terms  of  a  just 
detectable  target.   (A  precise  definition  of  a  just  detectable 
target  can  be  made  in  terms  of  a  specified  detection  probability, 
false  alarm  probability  and  target  dipole  moment.)   The  default 
choice  for  MAXF  is  consistent  with  the  observation  in  Reference  1 
that  if  an  optimum  value  for  MAXF  is  determined  for  a  minimum 
dipole  moment  target,  then  no  significant  increase  in  MAXF  is 
required  in  order  to  maintain  a  required  detection  probability  if 
the  encounter  lateral  range  is  decreased  even  though  the  signal 
energy  spectrum  is  shifted  to  higher  frequencies. 

Because  of  the  detection  statistics  and  the  gaussian  noise 
model,  if  there  were  no  penalty  for  decision  delay,  a  sample 
interval  length  for  a  signal  should  be  chosen  equal  to  the  signal 
duration,  since  this  would  make  the  detection  probability  for  an 
encounter  a  maximum.   In  the  program,  a  default  choice  for  the 
sample  interval  length  (integration  time)  is  2 • MAXRO/MINVM  where 
MINVM  is  a  user  estimated  maximum  encounter  relative  speed 
converted  from  knots  to  meters  per  second  and  MAXRO  is  a  user 
estimated  maximum  slant  range  at  CPA  in  meters  for  a  detectable 
encounter  (in  terms  of  a  specific  detection  probability  and  false 
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alarm  probability)  with  a  user  specified  maximum  dipole  moment 
target.   This  choice  might  be  considered  a  balance  between 
minimizing  decision  delay  and  maximizing  detection  probability. 
In  the  encounter  models,  the  sample  intervals  are  located  so  that 
the  CPA  time  is  at  the  center  of  a  sample  interval  and  this 
sample  interval  is  the  only  one  that  contains  signal  energy. 
This  characteristic  is  consistent  with  the  default  choice  for  the 
sample  interval.   For  a  given  sample  interval  length,  although  in 
general  the  interval  location  is  not  the  optimum  one  in  terms  of 
signal  energy,  it  should  be  approximately  so  in  most  cases. 
The  program  noise  parameter  is  SIG.   It  represents  the 
standard  deviation  a      associated  with  the  magnetic  noise  process 
of  the  two  encounter  models.   In  terms  of  the  ideal  low  pass 
filter  implied  by  the  encounter  models,  its  sguare  should  be 
egual  to  MAXF • (SIGO) 2  where  SIGO  is  the  magnetic  noise  process 
constant  power  spectral  density.   The  program  does  not  enforce 
this  relation.   Therefore,  in  using  the  program,  the  implied 
relation  between  the  two  input  parameters:  magnetic  noise  and 
maximum  magnetic  signal  freguency  should  be  kept  in  mind.   If  an 
average  value  of  the  peak-to-peak  magnetic  noise  for  an  encounter 
can  be  estimated,  for  example  from  a  magnetometer  trace,  then  the 
value  for  SIG  should  be  chosen  so  that  the  estimate  is  4  to  6 
times  this  value. 
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Appendix  1.   The  Detection  Statistics 

In  this  appendix,  y1#  y2 ,  •••  ,    Ym      are  sequential  values 
(voltages)  representing  the  sample  values  in  a  sample  interval. 
They  are  the  input  to  a  magnetometer's  detector.   With  these  m 
sample  values,  the  detector  computes  the  value  of  a  detection 
statistic.   This  value  is  represented  by   x   and  the  detection 
decision  corresponding  to  the  sample  interval  is  determined  by 
the  decision  rule:  If  x  >  x* ,  then  the  input  during  the  sample 
interval  was  noise  plus  signal,  otherwise,  the  input  was  noise. 

For  both  encounter  models,  the  detection  probability  and  the 
false  alarm  probability  are  decreasing  functions  of   x    and  the 
relation  is  one-to-one  in  both  cases.   In  the  program,  the  false 
alarm  probability  Pf   is  used  to  determine  a  unique  value  of 
the  threshold   x* .   This  value  is  then  used  to  determine  a  unique 
value  of  the  detection  probability  p^. 

In  the  program,   pf   is  found  using  the  relation  pf  =  R-St 
where   R   is  the  false  alarm  rate  in  false  alarms  per  second 
and   &t   is  the  sample  interval  length  in  seconds.   This  relation 
is  based  on  the  following  argument:   With  no  signal  energy  in  a 

sample  interval,  y±   =  n-±,    Y2   ~  n2'     '  '  '  »  Vm  =  nm  where   n^,  n2 , 
• • •  ,  nm   are  noise  values  (voltages)  input  to  a  magnetometer's 
detector.   The   n^   are  values  of  independent  normal  (gaussian) 
random  variables,  each  with  mean  zero  and  standard  deviation   a. 
Because  of  this,  in  the  encounter  models,  values  of   x   for 
different  sample  intervals  are  the  values  of  independent  random 
variables  that  determine  two  outcomes:   x  >  x*   or   x  <  x*. 
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Therefore,  in  terms  of  these  outcomes,  an  encounter  is  a  sequence 
of  independent  Bernoulli  trials.   If  there  is  no  signal  present, 
since  the  noise  is  determined  by  a  stationary  process,  pf   is  the 
same  for  each  sample  interval  and,  in  terms  of  these  outcomes, 
the  sequence  is  a  series  of  repeated  independent  Bernoulli  trials 
and  therefore  the  expected  number  of  trials  between  false  alarms 
is   1/Pf.   Since   the  time  between  trials  is   St,  the  expected 
number  of  seconds  between  false  alarm  is   6t/pf   or  the  expected 
number  of  false  alarms  per  second  R   is   equal  to   pf/St. 

The  determination  of   x*   depends  on  the  encounter  model 
statistic.   For  both  encounter  models,  when  there  is  a  signal, 
Yl  =  nl  +  slf  y2  =  n2  +  y2,  ■ • •  ,  ym  =  nm  +  sm  where   slf  s2, 
•••  ,  sm   are  signal  values  (voltages)  input  to  a  magnetometer's 
detector.   The  models  imply  that  the  signal  values   s-[  =  K-(HS)^ 
where   K   is  a  constant  whose  value  is  determined  by  the 
characteristics  of  the  encounter  magnetometer  and  where  the 
(Hs)^   are  dipole  magnetic  signal  intensities.   The  models  also 
imply  that  the  noise  values   n^   are  determined  by  a  gaussian 
stochastic  processes  characterized  by  a  standard  deviation   a 
and  that  n-[  =  K-  (HN)  j_  +  n^  where  the   (%j)  i   are  magnetic  noise 
intensities  and  the   n-^   are  magnetometer  instrument  noise  values 
(voltages).   In  the  program,  the  magnitude  of   K   is   1.   Since, 
for  both  models,   p^   depends  only  on  the  ratio  of  signal  energy 
to  noise  energy  for  a  sample  interval,  this  is  a  satisfactory 
choice  for  the  program  if  instrument  noise  is  assumed  to  be 
determined  by  a  process  independent  of  the  magnetic  noise  process 
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and  to  be  expressed  in  terms  of  an  equivalent  magnetic  noise 
(H)l  =  (1/K)  n[. 

For  both  encounter  models,  the  signal  (the  average  signal 
power)   S  =  (1/m)  S  s^   and  the  noise  (the  expected  value  of  the 
average  noise  power)   N  =  a    so  that  the  signal-to-noise  ratio 
is   (1/m) -2  s^/a2   where  the  sum  index   i  =  1,  2,  • • •  ,  m. 
The  Crosscorrelation  Detector  Statistic:   The  statistic  for  the 
first  encounter  model  is  a  crosscorrelation  detector  statistic 
that  is  defined  by  the  sum 

x  =  Z  Yi-Si 
where  the  summation  index   i  =  1,  2,  ■ • ■     ,  m   and  the  sum  is  over 
the  values  corresponding  to  a  sample  interval.   For  the  first 
encounter  model,  the  characteristics  of  both  the  noise  and  the 
signal  are  required  in  order  to  determine  encounter  detection 
probabilities.   In  particular,  the  signal  values  for  an  encounter 
are  in  the  memory  of  the  detector  prior  to  the  encounter.   For 
the  encounter  conditions  and  a  specified  false  alarm 
probability,  the  statistic  is  optimum  in  the  sense  that  the 
encounter  detection  probability  for  this  statistic  is  at  least 
equal  to  that  for  any  other  statistic.   Because  of  these 
considerations,  encounter  probabilities  based  on  the 
crosscorrelation  statistic  can  be  considered  to  represent  upper 
bounds  on  detection  performance  against  dipole  targets  for 
magnetometers  of  the  type  described  by  the  models. 

For  a  sample  interval  without  signal  energy,  x   is  the  value 
of  a  normal  random  variable  with  a  mean   (ix  =  0   and  a  variance 
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&X  =   °    '2  si   where   a   is  the  standard  deviation  associated  with 
the  noise  process  and  the  sum  index   1=1,2,  • • •  ,  m   and  the 
sum  is  over  the  values  corresponding  to  the  sample  interval. 
This  implies  that 

pf  =  1  -  P(x*/ax) 
where   P(z)   is  the  standard  normal  cumulative  distribution 
function.   This  relation  is  the  basis  for  determining  the 
threshold  value   x*. 

For  the  sample  interval  with  signal  energy,  x   is  the  value 

2 

of  a  normal  random  variable  with  a  mean   n-x  =  2  si   where  the  sum 
index   i  =  1,  2,  • ••  ,  m   and  the  sum  is  over  the  values 
corresponding  to  the  sample  interval.   This  implies  that 

pd  =  1  -  P(v*  -  d*) 
where   v*  =  x*/ax   and   d  =  S  s-j/a2  =  (l/m)-S/N.   This  relation 
is  the  basis  for  determining  encounter  detection  probabilities 
for  the  first  encounter  model.   The  relation  implies  that  for  a 
specified  false  alarm  probability  pf   the  detection  probability 
P3  is  an  increasing  function  of  the  signal  to-noise  ratio   S/N. 
The  Square  Law  Detector  Statistic:   The  statistic  for  the  second 
encounter  model  is  a  square  law  (energy)  detector  statistic  that 
is  defined  by  the  sum 

x  =  Z  y[ 
where  the  sum  index   i  =  1,  2,  • • •  ,  m   and  the  sum  is  over  the 
values  corresponding  to  the  sample  interval.   For  the  second 
encounter  model,  only  the  characteristics  of  the  noise  are 
required  to  determine  encounter  detection  probabilities. 
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For  a  sample  interval  without  signal  energy,  x/a2   is  the 
value  of  a  chi-square  random  variable  with  m  degrees  of 
freedom.   This  implies  that 

pf  =  1  -  P(x*/o2\m) 
where   P(x*/a2|m)   is  the  chi-square  cumulative  distribution 
function  for  a  chi-square  random  variable  with  m  degrees  of 
freedom  and  where  a      is  the  standard  deviation  associated  with 
the  noise  process.   This  relation  is  the  basis  for  determining 
the  threshold  value   x  . 

2 

For  the  sample  interval  with  signal  energy,  x/a    is  the 
value  of  a  noncentral  chi-square  random  variable  with  m 

2      2 

degrees  of  freedom  and  noncentral  parameter   £  s^/a   where  the 
sum  index   i  =  1,  2,  • • •  ,  m   and  the  sum  is  over  the  values 
corresponding  to  the  sample  interval.   This  implies  that 

pd  =  1  -  P(x*/°2|m,  S  Si/a2) 
where   P(x  /a  |m,  Z  s^/a    )   is  the  noncentral  chi-square 
cumulative  distribution  function  for  a  noncentral  chi-square 
random  variable  with  m  degrees  of  freedom  and  noncentral 

2      2 

parameter   E  s^/a  =  (l/m)-S/N.   This  relation  is  the  basis  for 
determining  encounter  detection  probabilities  for  the  second 
encounter  model.   The  relation  implies  that  for  a  specified  false 
alarm  probability  pf,  the  detection  probability  pd  is  an 
increasing  function  of  the  signal-to-noise  ratio   S/N.   This  is 
made  more  evident  by  the  following  relation: 

P(x*/</|m,  s  s^/a2)  =  S{  (aJ/J  ! )  ■  exp(-a)  •  P[x*/°2  I  (m  +  2-j)]} 
where  the  sum  index   j  =  0,  1,2,  • • •   and  the  parameter 
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a  =  (1/2)  Z  Sj/a  =  (m/2) • (S/N)  where  the  sum  index  1=1,  2, 
• • •  ,  m.  Note  that  the  relation  p^  >  pf  is  satisfied,  since 
P(x*/a2|m)  <  P[x*/a2|(m  +  2j)]   for   j  =  0,  1,  2,  •••  . 
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Appendix  2.   Program  Probability  Calculations 

The  program  evaluates  the  cumulative  and  inverse  cumulative 
distribution  functions  using  approximations  described  in 
Reference  2.   These  approximations  are  listed  below. 
The  Standard  Normal  Cumulative  Distribution  Function 
Approximation : 

P(z)  =  1  -  s-t-(b!  +  t- (b2  +  t-(b3  +  t- (b4  +  t-b5)))) 
where   s  =  (1/  (2  •  tt)  *)  •  exp(-z  2/2)      and  t  =  1/(1  +  b0-z).   And 
where 

b0  =  .2316419,   bx  =  .319381530,   b2  =  -.356563782, 

b4  =  -1.821255978,   b5  =  1.330274429, 
and   z  >  0.   For   z  <  0,  P(z)  =  1  -  P(|z|). 

The  Inverse  Standard  Normal  Cumulative  Distribution  Function 
Approximation : 

z(P)  =  t  -  (c0  +  t(Cl  +  tc2))/(l  +  t(dx  +  t(d2  +  td3))) 
where   t  =  (ln(l/Q2))^   and   Q  =  1  -  P.   And  where 

c0  =  2.515517,  c±    =    .802853,   c2  =  .010328, 

dx  =  1.432788,   d2  =  .189269,   d3  =  .001308, 
and   .5  <  P  <  1.   For   0  <  P  <  .5,   z(P)  =  -z(l  -  P) . 
The  Inverse  Chi-Square  Cumulative  Distribution  Function 
Approximation : 

v(P|m)  =  m-[l  -  2/(9-m)  +  z • (2/ (9 -m) ) \ ] 3 
where   P(z)  =  P(v|m).   In  the  program,  the  inverse  standard 
normal  cumulative  distribution  function  approximation  is  used  to 
determine   z. 
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The  Noncentral  Chi-Square  Cumulative  Distribution  Function 
Approximation : 

P(w|m,Z  Sj/a  )  =  P(z) 
where   z  =   [2-w/(l  +  b)]i   -  [2 -a/ (1  +  b)  -  1]  §   with 

2      2  2      2  2      2 

a  =  m  +  E  s^/a  ,   b  =  (Z  s-jya  )/(m  +  2  s^/a  )   and  the  sum  index 
i  =  1,  2,  • • •  ,  m.   In  the  program,  the  standard  normal 
cumulative  distribution  function  approximation  is  used  to 
determine   P(z) . 
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Appendix  3.   The  Magnetic  Signal 

The  encounter  models  are  defined  by  the  following 
conditions:   A  submarine  magnetic  anomaly  field  is  a  magnetic 
dipole  field  that  is  superimposed  on  an  earth  magnetic  field  that 
is  constant  over  an  encounter  region.   A  magnetometer  magnetic 
signal  value  is  the  magnitude  of  the  projection  of  a  dipole 
magnetic  field  on  the  earth  magnetic  field  at  the  location  of  the 
magnetometer.   The  basis  for  determining  magnetic  signal  values 
in  the  program  is  an  expression  that  can  be  developed  as  follows: 
In  the  right  handed  rectangular  coordinate  system  that  is  shown 
in  Figure  1,  a  magnetic  dipole  is  at  the  origin,  the  xy-plane  is 
the  horizontal  plane  at  a  representative  point  in  an  encounter 
region,  the  positive  y-axis  is  in  the  direction  of  magnetic 
north,  the  positive  x-axis  is  in  the  direction  of  magnetic  east 
and  the  positive  z-axis  is  positive  upward.   In  this  rectangular 
coordinate  system,  the  constant  earth  magnetic  field  can  be 
expressed  by   He  =  He(jcos  $  +  k-sin  §)   where   $   is  a 
magnetic  dip  angle  and  He   is  a  magnetic  field  magnitude  that 
characterizes  the  earth  field  in  an  encounter  region. 

In  a  spherical  coordinate  system  with  the  origin  at  the 
magnetic  dipole  and  the  polar  axis  in  the  direction  of  the  dipole 
moment,   Hd  =  (c -p/r3) • (2 • r- cos  6  +  e-sin  6)   is  the  magnetic 
field  of  the  dipole  at  a  point  with  spherical  coordinates 
(r,r,9) .  In  this  expression   p   is  the  magnitude  of  the  dipole 
moment  and   c   is  a  constant  whose  value  is  determined  by  the 
choice  of  units.   The  magnetic  signal  for  a  magnetometer  that  is 
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Figure  1.  A  unit  vector  in  the  direction  of  the  earth  field  and  a 
unit  vector  in  the  direction  of  the  dipole  field  of  a 
dipole  at  the  origin  are  shown  at  a  point  that  is  a 
distance   r   from  the  dipole.   The  unit  vector  h   is 
in  the  direction  of  the  horizontal  component  of  the 
dipole  moment. 
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described  by  the  encounter  models  is   Hs  =  He'H^/He  when  the 
magnetometer  is  at  the  point  (r,r,0)  and  Hs   can  be  expressed  in 
terms  of  the  rectangular  coordinate  system  as  follows:   First, 
let   r0  be  a  unit  vector  in  the  direction  of  the  spherical 
coordinate  system  polar  axis.   Then  the  magnetic  dipole  field 

Hd  =  (c-p/r3) • (3-r-cos  6  -  r0) 
since   r  =  (ix  +  jy  +  kz)/r  where   r  =  (x*  +  y2  +  z2)1    , 
6  =  [ (r0  x  r)  x  r]/sin  G  =  [(r-r0)r  -  (r-r) • r0]/sin  6,   and  the 
dot  product   r»r0  =  cos  6.   The  unit  vector   r0   can  be  expressed 
in  the  rectangular  coordinates  by  noting  that   r0  =  p/p   and  then 
expressing  p   in  rectangular  coordinates.   To  do  this,  let  n 
be  the  depression  angle  of   p   from  the  xy-plane  (the  horizontal 
plane)  with   n   positive  downward  and  let   a   be  the  direction  of 
p  relative  to  magnetic  north.   Then  the  magnetic  dipole  moment 
p  =  p-  (h-cos  n  -   k-sin  n)   in  terms  of  the  unit  vector   k  and 
the  unit  vector   h  =  i-sin  a  +  j • cos  n  which  has  the  direction 
of  the  horizontal  component  of  p.   With  these  results,  the  unit 
vector   r0  =  i-  (cos  n-sin  a)  +  j  •  (cos  n-cos  a)  -  k-sin  n   and 
Hs  =  (c-p/r3) •[ (3/r) -cos  9(ycos  $  -  z • sin  $) 

-  (cos  $  •  cos  n-cos  a  +  sin  $-sin  n)  ] 
where   cos  9  =  (l/r)-(xcos  n-sin  a  +  y-cos  n-cos  a  -  z • sin  n) 
since   cos  9  =  r-r0.   As  can  be  seen  from  this  expression,  for  a 
constant  dipole  moment  magnitude  and  direction  and  a  constant 
earth  field  magnitude  and  direction,  the  magnetic  signal  is  only 
a  function  of  the  rectangular  coordinates  of  the  location  of  the 
magnetometer  relative  to  the  dipole.   In  the  encounter  models, 
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both  of  these  conditions  are  satisfied.  However,  by  allowing  p, 
n,  a,  Hs  and  $  to  vary,  the  expression  for  Hs  is  applicable 
to  more  general  encounter  models. 
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Appendix  4 .   The  Anderson  Formulation 

In  the  encounter  models,  the  magnetic  signal   Hs   at  a 
sample  point  in  a  straight  line  encounter  can  be  represented  in  a 
form  described  by  Anderson  in  Reference  4.   For  convenience,  the 
Anderson  formulation  is  used  in  the  program  to  determine  values 
for  Hs.   It  can  be  developed  as  follows:   In  Figure  2,  the 
primed  rectangular  coordinate  system  is  superimposed  on  the 
rectangular  coordinate  system  of  Figure  1  so  that  the  origin  is 
coincident  with  the  origin  of  that  system.   A  magnetometer  is  in 
a  straight  line  encounter  with  a  magnetic  dipole  that  is  located 
at  the  origin  of  the  combined  system.   The  combined  system  moves 
with  the  magnetic  dipole  with  the  x'-axis  oriented  so  that  it  is 
parallel  to  and  in  the  direction  <p      of  magnetometer's  relative 
to  the  magnetic  dipole  and  the  z'-axis  oriented  so  that  it  is 
directed  toward  and  passes  through  the  CPA  on  that  track.   Let 
1,  m,  and   n   be  the  direction  cosines  of  the  dipole  moment   p 
and   1]_,  n^   and   n^   be  the  direction  cosines  of  the  earth 
magnetic  field  He  in  the  primed  system.   Then  the  unit  vector 
Tq  =  i'-l  +  j'n  +  k' -m.   And,  for  points  on  the  relative  track, 
x'  =  s1,  y*  =  0   and   z1  =  R  where   R   is  the  slant  range  of  the 
dipole  at  CPA   and  where   s'  is  the  algebraic  distance  of  the 
magnetometer  from  CPA  on  the  relative  track.  (It  is  negative  for 
points  before  CPA  and  positive  for  points  after  CPA.)   This 
implies  that  magnetometer  position  vector  in  the  moving 
coordinate  system  is   r  =  i1  (s'/r)  +  k' ■ (R/r) . 
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Figure  2 .  The  dipole  is  at  the  origin  in  the  unprimed  and  primed 
coordinate  systems.   In  the  primed  coordinate  system, 
the  CPA  is  at   (0,  0,  R)   and,  for  a  time   t   during 
the  encounter,  the  magnetometer  is  at   [s* (t) ,  0,  R] . 
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From  Appendix  3 , 

Hd  =  (c-p/r3) ■ (3 -r- cos  6  -  r0) • 
With   r  and   r0   expressed  in  terms  of  the  primed  unit  vectors 
and   cos  0  =  r«r0  =  l(s'/r)  +  n(R/r),  this  becomes 

Hd  =  (c-p/r3) •[ (3/r2) • (Is'  +  n-R)-(i'-s'  +  k'  -R) 
-  (i'  1  +j'  m  +  k' -n) ]. 
Then,  since    He/He  =  ±,11   +   j'-n^  +   k'-nx   and   Hs  =  He-Hd/He, 
Hs  =  (cp/r3) • [ (2-l-lx  -  m-mx  -  nn1)(s')2 

+  3-(nlx  +  l-n^J-s'-R  +  (2nnx  -  l-li  -  m-mxJ-R2] 
The  quantities 

A2  =  2'1'ln  -  m-m^  ~  n-n^ 

A±    =    3  •  (n-  1-l  +  1  -n^ 

A0  =  2-n-n^  -  lli  -  m-m^, 

1 
are  the  Anderson  coefficients.   With   r  =   [(s')   +  R  ]2   and 

P  =  s'/R,   Hs   can  now  be  expressed  as  follows: 

HS(P)  =  (cp/R3) -2  Aj-Fj(P) 

where  the   Fj(p)   =   pj/(l  +  P2)5//2   for   j  =  0,  1,  2   are  the 

Anderson  functions.   This  is  the  Anderson  formulation. 

To  relate  the  Anderson  formulation  for  Hs   to  the 

formulation  for  Hs   in  Appendix  3,  first  note  that 


1 

= 

ro-i' 

m 

= 

r0-j' 

n 

= 

r0-k' 

and 
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1-L  =  (He/He)  -i' 
m±   =    (He/He) -j' 

nx  =  (He/He) -k' . 

Then  express   r0/  (He/He)   and  the  unit  vectors   i' ,  j*  and   k' 
in  terms  of  the  unit  vectors   i,  j   and   k  and  take  the 
indicated  dot  products.   From  Appendix  3, 

r0  =  i (cos  fi-sin  a)  +  j- (cos  ft • cos  a)  -  k-sin  ft 
and 

He/He  =  j • cos  $  +  k-sin  $. 
To  express   i1 ,  j'  and   k'  in  terms  of  the  unit  vectors   i,  j 
and   k,  note  that  the  unprimed  coordinate  system  can  be 
transformed  to  the  primed  coordinate  system  by  two  rotations  that 
are  defined  as  follows:   First,  rotate  a  coordinate  system  that 
is  coincident  with  the  unprimed  coordinate  system  about  the 
z-axis  through  the  angle  (0  -  ir/2)  with  positive  angles  clockwise 
(left  hand  rule)  so  that  its  x-axis  is  parallel  to  and  in  the 
direction  of  the  relative  track.   Next,  rotate  this  system  about 
its  x-axis  through  an  angle   5   with  positive  angles 
counterclockwise  (right  hand  rule)  so  that  the  positive  z-axis 
passes  through  the  CPA.   The  angle   6   is  related  to  the 
vertical  separation   z0   of  the  magnetometer  and  the  dipole  and 
the  algebraic  encounter  lateral  range   L  that  is  positive  if  the 
dipole  is  to  the  left  of  the  relative  track.   With  these  sign 
definitions:   L  =  R-sin  S   and   z0  =  R-cos  6.   After  the 
rotation,  the  auxiliary  coordinate  system  is  coincident  with  the 
primed  coordinate  system. 

27 


These  transformations  can  be  described  in  terms  of  matrix 
equations  as  follows:   Let   (x",y",zM)   be  the  coordinates  of  a 
point  in  the  coordinate  system  that  is  coincident  with  the 
auxiliary  coordinate  system  after  the  first  rotation.   Then  the 
transformation  from  the  unprimed  coordinate  system  to  this  double 
primed  coordinate  system  is  described  by  the  matrix  equation 


X" 

y" 

= 

z" 

sin  0    cos  0    0 

•cos  0    sin  0    0 

0        0        1 


x 

y 

z 


And  the  transformation  from  the  double  primed  coordinate  system 
to  the  primed  coordinate  system  is  described  by  the  matrix 
equation 


X' 

y' 

z' 


X" 

y" 

Z" 

10         0 

0      cos  5     sin  5 

0     -sin  S     cos  & 

Taking  the  product  of  the  rotation  matrices  in  the  indicated 
order  yields  the  matrix  equation 

sin  0  cos  0  Ox 

■cos  5  cos  0    cos  5  sin  0    sin  6     y 
sin  5  cos  0    -sin  5  sin  <p  cos  5      z 

which  defines  the  transformation  from  the  unprimed  to  the  primed 
coordinate  system.   The  unprimed  vector  components  of  the  unit 
vector   i»  can  be  found  by  transforming  the  coordinates  (1,0,0) 


X1 

y 

= 

z' 
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in  the  primed  system  to  their  corresponding  coordinates  in  the 
unprimed  system  with  the  inverse  of  this  matrix  and  then 
repeating  this  process  for  (0,1,0)  and  (0,0,1)  in  order  to  find 
the  unprimed  unit  vector  components  of   j1   and   k' .   However, 
since  the  inverse  transformation  matrix  is  the  transpose  of  this 
matrix,  the  elements  of  the  row  that  corresponds  to  a  primed  unit 
vector  are  the  magnitudes  of  the  unprimed  vectors  that  are  its 
components.   Consequently: 

i1  =   i-sin  0  +  j • cos  0 

j •  =  -i-cos  S  cos  0  +  j • cos  8  sin  0  +  k-sin  & 
k'  =   i-sin  6  cos  0  -  j • sin  5  sin  0  +  k-cos  S. 
Then,  taking  the  dot  products  between   r,  r0   and  these  three 
unit  vectors  as  indicate  above  gives: 

1  =  cos  ficos  (0  -  a) 

m  =  cos  5  •  cos  fi-sin  (0  -  a)  -  sin  S-sin  f2 
n  =  -sin  5 • cos  n-sin  (0  -  a)  -  cos  6-sin  n 
and 

1^  =  cos  § ■ cos  0 

m^  =  cos  6 • cos  $ • sin  0  -  sin  & • sin  § 

n^  =  -sin  & • cos  fsin  0  -  cos  6-sin  §. 

These  are  the  relations  that  are  used  in  the  program  to  determine 
values  for  the  Anderson  coefficients. 
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Appendix  5.   The  Encounter  Equations  of  Motion 

In  the  double  primed  coordinate  system  that  is  defined  in 
Appendix  4,  the  equations  of  motion  of  a  magnetometer  relative  to 
a  submarine  (dipole)  are: 

x"  (t)  =  s'  (t) 

y--(t)  =  -L 

z»(t)  =  z0. 
where   L   is  the  algebraic  encounter  lateral  range  that  is 
defined  in  Appendix  4,  z0   is  the  vertical  separation  between  the 
magnetometer  and  the  submarine  and   s'(t)   is  the  distance  of  the 
magnetometer  from  the  CPA  on  the  relative  track.   With  w   the 
speed  of  the  magnetometer  relative  to  the  submarine  and   t   a 
relative  time  parameter,   s'(t)  =  wt  .   These  equations  can  be 
considered  to  be  the  ones  used  in  the  program  to  describe  the 
motion  of  a  magnetometer  relative  to  a  submarine.   There,   t   is 
determined  by   t  =  [j  -  (m-l)/2]-5t  where  the  index   j  =  1,  2, 
• • •,  m   and   5t,  a  time  step,  is  the  time  between  samples.   Note, 
when   t  =  0,  the  magnetometer  is  at  the  CPA. 

In  the  coordinate  system  of  Figure  1  in  a  straight  line 
encounter  as  defined  in  the  encounter  models,  the  equations  of 
motion  of  a  magnetometer  relative  to  a  submarine  can  be  written 
as  follows: 

x(t)  =  s'(t)sin  0   +   L-cos  0 

y(t)  =  s'(t) -cos  0   -   L-sin  0 

z(t)  =  z0, 
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since  the  transformation  from  the  double  primed  coordinate  system 
to  the  primed  coordinate  system  is  determined  by  the  matrix 
equation 


x 

y 

z 


sin  0    -cos  0    0 
cos  <p  sin  0    0 

0        0       1 


X" 

y" 

z" 

x"  =  s1 (t)   and  y"  =  -L.   The  above  equations  and  the  expression 
for   Hs   in  Appendix  3  could  have  been  used  in  the  program  to 
evaluate  the  magnetic  signal.   In  particular,  with  these 
equations  of  motion  and  with  the  two  relations   L  =  R-sin  S      and 
Zq  =  R-cos  6,  the  expression  for  Hs   in  Appendix  3  can  be 
written  in  terms  of   ft,  a,  $,  0,  5   and   R   so  that  it  is 
identical  in  appearance  to  the  Anderson  formulation  for   Hs   in 
terms  of  these  quantities.   The  definition  of   &   in  Appendix  4 
in  terms  of  a  counterclockwise  rotation  results  in  a  definition 
of  the  algebraic  lateral  range  that  is  consistent  with  some  that 
have  been  used  elsewhere. 

In  the  program,  the  relation  w  =  v  -  u   is  the  basis  for 
determining  the  relative  speed   w.   In  this  relation,  v   is  the 
velocity  of  the  magnetometer,  u   is  the  velocity  of  the  submarine 
(dipole)  and  w   is  the  velocity  of  the  magnetometer  relative  to 
the  submarine.   This  relation  implies  the  following  equations: 
wx  =  vx  -  ux,  Wy  =  Vy  -  Uy   and  wz  =  vz  -  uz   where  the 
coordinates   x,  y   and   z   refer  to  a  fixed  coordinate  system 
with  the  same  orientation  as  that  of  Figure  1.   In  the  encounters 
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of  the  models,  vx  =  v-sin  a,  vy  =  v-cos  a   and  vz  =  0   where  a 
is  the  magnetic  course  and  v   is   speed  of  the  magnetometer. 
And,  ux  =  u-sin  (3 ,  uy  =  u-cos  (3   and  uz  =  0   where   (3   is  the 
magnetic  course  and  u   is  the  speed  of  the  submarine.   The 
relative  magnetic  course  0   and  the  relative  speed  of  the 
magnetometer  are  defined  by   wx  =  w-sin  0   and   wy  =  w-cos  0.   In 
the  program,  <p      and  w  are  determined  with  a  rectangular  to 
polar  conversion  routine  where  w  =  (  wx  +  wy) 5   and  where   0   is 
determined  by   sin-1(wx/w)   and   cos-1(wy/w). 
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Appendix  6.   The  Submarine  Magnetic  Dipole 

In  the  encounter  models,  the  magnitude  and  direction  of  a 
submarine's  dipole  moment  p   are  determined  by  first  determining 
its  components  in  the  rectangular  coordinate  system  of  Figure  1. 
In  that  coordinate  system, 

Px  ~  (PLP  +  PLl)  sin  P  +  (PTP  +  PTl)  cos  P 

Py  =  (PLP  +  PLl) 'cos  P  "  (PTP  +  PTl) 'sin  P 

Pz  =  ~(PVP  +  Pvi) 
where   (3   is  the  submarine's  magnetic  course,  Plp,  Ptp   ^rid  pyp 

are  the  permanent  and   PlI/  Pti   an<^  Pvi   are  the  induced 
longitudinal,  transverse  and  vertical  magnetic  dipole  moments  of 
the  submarine.   These  relations  are  based  on  the  following  sign 
convention: 

pL  is  positive  when   pL   is  directed  from  stern  to  bow. 

pT  is  positive  when   pT   is  directed  from  port  to  starboard. 

Py  is  positive  when  py   is  directed  downward. 
The  permanent  dipole  moments  are  input  parameters  in  the  program 
and  the  induced  dipole  moments  are  determined  in  a  way  that  is 
similar  to  one  that  is  described  in  Reference  1.   In  the 
encounter  models,  a  submarine  is  a  ferromagnetic  prolate 
ellipsoid  with  the  major  axis  the  submarine's  longitudinal  axis 
and  the  equal  minor  axes  the  submarine's  transverse  and  vertical 
axes.   And,  the  induced  dipole  moments  are: 

PLI  =  kLHeL 
PTI  =  kT " HeT 

Pvi  =  kvHev 
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where   HeL,  HeT   and   HeV  are  the  vector  components  of  the  earth 
magnetic  field  in  the  rectangular  coordinate  system  defined  by 
the  ellipsoid  axes  and  the  submarine  magnetic  moment  sign 
convention.   The  magnitude  of  these  vector  components  are: 

HeL  =   Hecos  §  cos  p 

HeT  =  -He  •  cos  $  sin  (3 

HeV  =  He' sin  $ • 
The  earth  magnetic  field  dip  angle   §   and  the  earth  magnetic 

field  magnitude   He   can  each  be  chosen  to  be  input  parameters  in 

the  program  or  the  can  be  computed  by  the  program  as  described  in 

Appendix  7 .  (Values  for   $   and  He   can  be  found  from  magnetic 

data  charts,  for  example,  see   References  5  and  References  6.) 

By  using  the  above  relations, 

Px  =  He-cos  §•  (kL  -  kT)  -sin  (3  cos  p  +  (pLP-sin  P  +  pTP-cos  P) 

Py  =  He-cos  $(kLcos2p  +  kT-sin2p)  +  (pLPcos  p  -  pTP-sin  P) 

pz  =  -(Hekv-sin  $  +  pvp) 

In  the  encounter  models,  the  values  of  the  permeability 

coefficients   kL,  kT  and  kv  are  related  to  submarine 

displacement  by  the  following  relations: 

kL  =  fL-W 

kT  =  f T • W 

kv  =  fv-w 
where   W   is  the  submarine  displacement  in  tons  and   fL,  fT   and 
fv  are  permeability  factors  that  are  determined  by  a  submarine's 
magnetic  characteristics.   In  the  program,  the  units  of  He   are 
gamma,  the  units  of   p   are  in  oersted-centimeter3,  the  units  of 
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k  are  oersted-centimeter3/gamma  and  the  units  of   f   are 
oersted-centimeter3/gamma-ton.   If  the  units  of   He   were 
gamma,  but  the  units  of   p  were  gamma-foot3,  then  the  units  of 
k  would  be  foot3  and  the  units  of   f  would  be  foot3/ton.   To 
convert   p   in  gamma-foot3  to  oersted-centimeter3  or   k   in 
foot3  to  oersted-centimeter3/gamma  or   f   in  foot3/ton  to 
oersted-cent imeter3/gamma-ton,  divide  by  3.53.  (The  program 
default  values  for   fL,  fT   and   fy   are  values  from  Reference  1 
in  foot3/ton  that  have  been  divided  by  3.53  to  give  values  in 
oersted-centimeter3/gamma-ton. )   To  convert   p   in  weber-meter  to 
oersted-centimeter3,  multiply  by   [  1/ (4  •  tt)  ]  •  1010   and  to  convert 
p   in  ampere-centimeter2  to  oersted-centimeter3,  multiply  by 
1-103. 
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Appendix  7.   The  Earth  Magnetic  Field  and  Dip  Angle 

An  auxiliary  magnetic  field  model  is  described  in  this 
appendix.   The  model  is  the  basis  for  a  default  choice  for  either 
the  value  of  the  earth  magnetic  field  magnitude  parameter  He   or 
the  earth  magnetic  field  dip  angle  parameter   $.   Relative  to 
encounter  model  accuracy,  the  default  values  should  be  adequate 
in  most  cases. 

In  the  model,  the  earth  magnetic  field  is  generated  by  a 
magnetic  dipole  that  is  located  at  the  earth's  center  and  the 
earth  is  a  nonmagnetic  sphere  of  radius   re.   With  pe   the 
magnitude  of  the  dipole  moment,  the  magnitude  of  the  earth  field 
at  a  point  is 

He  =  (He-HeP  =  (cpe/r3)  •  (3cos26  +  1)  * 
where   9   is  the  polar  angle  of  the  point  in  a  spherical 
coordinate  system  and   c   is  a  constant  whose  value  is  determined 
by  the  choice  of  units.   The  dipole  moment  is  coincident  with  and 
in  the  direction  of  the  polar  axis  which  is  directed  toward  the 
earth's  southern  hemisphere.   In  this  coordinate  system,  at  any 
point  on  the  surface  of  the  earth: 

He  =  Heo-(3cos26  +  1)5 
where   Heo  is  the  value  of   He   at  the  magnetic  equator  which  is 
defined  by  the  points  on  the  earth  where  9  =  90°.   In  terms  of 
the  dip  angle,  at  any  point  on  the  earth's  surface: 

He  =  2-Heo- (3  cos2§  +  1)~*. 
This  expression  can  be  obtained  by  noting  that   $   can  be 
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defined  in  terms  of  the  or  the   r  and   6   components   of   He   in 
the  spherical  coordinates  as  follows: 

sin  $  =  -Her/He  =  -2-Heo- (cos  6)/He 
and 

cos  $  =  He0/He  =  Heo- (sin  G)/He. 
Based  on  these  relations, 

[(sin  $)/2]2  +  cos2$  =  (Heo) 2 
and 

He  =  Heo- [ (sin2)/4  +  cos's]"*  =  2 • Heo • (3 • cos 2$  +  1)~* 
The  dip  angle  is  determined  from  a  magnetic  latitude  for  the 
encounter  region. 

The  magnetic  latitude  and  longitude  of  the  representative 
point  of  an  encounter  region  can  be  defined  in  terms  of  its 
geographic  latitude  and  longitude  by  using  the  following 
transformations:   First,  convert  the  latitude  and  longitude  of 
the  point  to  rectangular  coordinates  in  a  right-handed  coordinate 
system  whose  origin  is  at  the  center  of  a  spherical  earth,  whose 
positive  z-axis  passes  through  its  north  geographic  pole  and 
whose  positive  y-axis  passes  through  the  point  on  its  eguator 
with  latitude  0°  and  longitude  0°.   Next,  rotate  a  coordinate 
system  that  is  coincident  with  this  system  in  a  clockwise 
direction  (left  hand  rule)  about  its  z-axis  so  that  its  positive 
y-axis  passes  through  the  point  with  latitude  0°  and  with 
longitude  equal  to  that  of  the  north  magnetic  pole.   Then,  rotate 
the  system  in  a  clockwise  direction  about  its  x-axis  so  that  its 
positive  z-axis  passes  through  the  north  magnetic  pole.   Next, 
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transform  the  rectangular  coordinates  of  the  representative  point 
in  this  system  to  its  coordinates  in  the  spherical  coordinate 
system  that  is  associated  with  it.   Then,  with  n  the  polar 
angle  of  the  representative  point  in  this  system,  the  magnetic 
latitude  of  the  representative  point  is   L^  =  90°  -  ft.   Since 
tan  §  =  -2  cot  6   from  above  and   6  =  180°  -   n,    the  dip  angle   $ 
is  given  by  the  following  relation: 

§  =  tan-1  (2  tan  Lj^)  . 
In  the  program,  the  transformations  described  above  are 
accomplished  in  part  by  a  rectangular  to  polar  conversion 
routine.   In  particular,  by  rotating  the  final  rectangular 
coordinate  system  about  its  z-axis  so  that  the  x-coordinate  is 
zero,  the  polar  angle  of  the  polar  coordinates  of  the 
representative  point  in  the  resulting  yz-plane  determines   Lj^. 

Some  values  of   $   and  He   that  are  listed  in  Table  1  were 
generated  using  the  program.   The  latitude  and  longitude  of  the 
magnetic  pole  that  are  in  the  program  are   76°  N  and  100°  W. 
Values  of   §   from  Reference  6  and  of   He   from  Reference  7  are 
also  listed  in  Table  1.   A  comparison  of  the  program  values  with 
these  values  gives  an  indication  of  the  errors  inherent  in  the 
procedure. 
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Encounter 

Encounter 

Program 

Chart  #3  0   Program 

Chart  #39 

Latitude 

Longitude 

Dip   $ 

Dip   $ 

Field  He 

Field  He 

60°N 

180°W 

74° 

70° 

.63 

.53 

60°N 

150°W 

78° 

73° 

.66 

.55 

60°N 

90°W 

82° 

84° 

.68 

.61 

60°N 

30°W 

75° 

74° 

.64 

.52 

30°N 

150°W 

58° 

50° 

.51 

.41 

30°N 

90°W 

62° 

60° 

.55 

.51 

30°N 

60°W 

59° 

59° 

.53 

.48 

30°N 

30°E 

37° 

43° 

.41 

.42 

30°N 

60°E 

31° 

45° 

.39 

.46 

30°N 

90°E 

30° 

44° 

.39 

.44 

30°N 

150°E 

42° 

40° 

.43 

.40 

0° 

90°W 

26° 

18° 

.38 

.32 

0° 

30°W 

9° 

-5° 

.35 

.28 

0° 

0° 

-48° 

-25° 

.35 

.31 

0° 

60°E 

-25° 

-20° 

.38 

.37 

30°S 

180°W 

-45° 

-55° 

.44 

.50 

30°S 

90°W 

-30° 

-34° 

.39 

.31 

30°S 

90°E 

-62° 

-66° 

.55 

.52 

60°S 

30°W 

-69° 

-56° 

.60 

.34 

Table  1.   Some  corresponding  values  for  the  encounter  dip  angle 
$   and  the  earth  magnetic  field  intensity   He   in 
gamma.   Chart  #33  refers  to  Reference  5  and  Chart  #39 
refers  to  Reference  6. 

A  small  difference  between  a  chart  value  and  a  computed  value  for 
the  dip  angle   $   at  some  location  does  not  imply  a  small 
difference  between  the  chart  value  and  the  computed  value  for  the 
earth  magnetic  field  intensity  at  that  location.   For  example, 
note  the  values  for   60°N   latitude  and   30°W   longitude. 
(Positive  values  of   $   indicate  the  inclination  or  dip  below  the 
horizontal  of  the  north  seeking  end  of  a  dip  needle.   Negative 
values  indicate  the  inclination  below  the  horizontal  of  the  south 
seeking  end. ) 
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An  extension  of  the  above  procedure  can  be  made  for  finding 
the  magnetic  variation  at  a  location.   However,  the  values 
generated  by  using  the  procedure  are  generally  unsatisfactory. 
Magnetic  variation  values  are  charted  in  Reference  7. 
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Appendix  8 .   An  Alternative  Encounter  Model 

The  alternative  model  that  is  referred  to  in  Section  III  is 
described  in  more  detail  in  this  appendix.   In  the  model,  the 
detection  range  is  the  slant  range   R   at  the  CPA  for  an 
encounter  with  a  specified  detection  probability  (usually  .5)  is 
defined  by: 

R  =  [c-p/Hs]V3 

where   c   and  p   are  defined  in  Appendix  3   and   Hs   represents 
a  minimum  detectable  average  magnetic  signal  that  is  defined  by: 

Hs  =  (ORF) -NM 
where   ORF  is  a  signal-to-noise  ratio  called  the  operator 
recognition  factor  and   NM   represents  the  magnetic  noise. 
Combining  these  two  relations  gives: 

R  =  <c-p/[(ORF)  -NM]}V3. 
The  value  for   ORF   depends  both  on  the  specified  detection 
probability  and  on  a  specified  or  implied  false  alarm 
probability. 

The  Anderson  formulation  is  consistent  with  these  relations 
in  an  approximate  sense  if  the  average  magnetic  signal   H   is 
defined  as  a  root  mean  sguare  value  such  that   H  =  (c-p/R3)  k 
where   k   is  an  encounter  parameter  defined  by: 


k  =  (E  [S  A-pF-jCPi)  ]2}" 


with  the  first  sum  index   i  =  1,  2,  • • •  ,  m   and  the  second  sum 
index   j  =  0,  1,  2.   For  a  particular  encounter  geometry,  k   is 
constant  and  this  suggests  that  the  two  encounter  models  could  be 
used  to  determine  an  average  value  for   k   for  an  encounter 
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region  based  on  average  submarine  magnetic  characteristics. 
Values  for  both  k  and  H  are  generated  by  the  program  and  such 
values  can  give  an  indication  of  the  magnitude  of  the  differences 
in  detection  range  estimates  that  are  based  on  this  model  and 
those  that  are  based  on  either  of  the  other  two  encounter  models. 
A  more  detailed  comparison  of  these  encounter  models  is  described 
in  Reference  8. 
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Appendix  9 .   An  Example  of  the  Program  Output 

The  program  is  designed  to  generate  the  following 
quantities:  encounter  parameter  values,  lateral  range  function 
values  for  the  crosscorrelation  encounter  model  and  for  the 
square  law  encounter  model,  average  magnetic  signal  values,  slant 
range  at  CPA  values,  encounter  parameter  values  and  magnetic 
signal  values.   These  values  can  be  saved  as  a  program  data  file 
and/or  they  can  be  printed. 

An  example  of  the  program's  printed  output  is  listed  in 
Table  2,  Table  3  and  Table  4.   Figure  3  is  a  plot  of  the  lateral 
range  function  values  (a  lateral  range  curve)  for  the  square  law 
detector  that  are  listed  in  Table  3.   Figure  4  is  a  plot  of  the 
magnetic  signal  values  that  are  listed  in  Table  4. 

A  second  example  of  the  program's  printed  output  is  listed 
in  Table  5,  Table  6  and  Table  7.   Figure  5  is  a  plot  of  the 
lateral  range  function  values  for  the  square  law  detector  that 
are  listed  in  Table  6.   Figure  6  is  a  plot  of  the  magnetic  signal 
values  that  are  listed  in  Table  7. 
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prograr,  file  name  MAD.BAS 

program  data  file  name  data. mad 

magnetic  data  file  name  data. mag 

processing  data  file  name  data. pro 

kinematic  data  file  name  data. kin 
combined  magnetic,  processing  &  kinematic  data  file  name  data.mpk 
encounter  latitude  (decimal  degrees)  45 

encounter  longitude  (decimal  degrees)  -45 

encounter  variation  (decimal  degrees)  -25 
encounter  dip  angle  (decimal  degrees)  68.41608 

encounter  magnetic  field  intensity  (gamma)  59035.05 

permanent  longitudinal  moment  (oersted-cm3)  0 

permanent  transverse  moment  (oersted-cm3 )  0 

permanent  vertical  moment  (oersted-cm3)  0 

target  displacement  (tons)  1800 

target  longitudinal  permeability  coefficient  13140 

target  transverse  permeability  coefficient  2880 

target  vertical  permeability  coefficient  2880 

target  longitudinal  permeability  factor  7.3 

target  transverse  permeability  factor  1.6 

target  vertical  permeability  factor  1.6 

sampling  period  (seconds)  .5 

integration  time  (seconds)  20 

adjusted  integration  time  (seconds)  20.5 

number  of  samples  per  encounter  41 

magnetometer  course  (decimal  degrees)  290 

magnetometer  speed  (knots)  180 

magnetometer  altitude  (meters)  100 

target  course  (decimal  degrees)  20 

target  speed  (knots)  10 

target  depth  (meters)  100 

magnetometer  relative  course  (decimal  degrees)  286.8202 

magnetometer  relative  speed  (knots)  180.2776 

magnetometer-target  vertical  separation  (meters)  200 

target  induced  longitudinal  dipole  moment  (oersted-cm3 )    2.017795E+08 
target  induced  transverse  dipole  moment  (oersted-cm3)     -4 . 422564E+07 
target  induced  vertical  dipole  moment  (oersted-cm3)        1.58099E+08 
magnetic  dipole  moment  (oersted-cm3)  2.601272E+08 

dipole  moment  azimuth  (decimal  degrees)  7.637512 

dipole  moment  depression  angle  (decimal  degrees)  37.42885 

distance  between  samples  on  the  relative  track  (meters)  46.3714 
false  alarm  rate  (false  alarms  per  hour)  3 

false  alarm  probability  1.666667E-02 

magnetic  noise  (gamma)  .35 

maximum  lateral  range  (meters)  500 

lateral  range  step  (meters)  20 

number  of  lateral  range  function  values  51 


Table  2.   An  example  of  an  encounter  parameter  values  printout 
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data.r.ad  lateral  range  function  values 


L 

p(cc) 

p(sl) 

H 

R 

k 

meters 

gamma 

meters 

-500 

.2155959 

.0255824 

.184247 

538. 

5165 

2.818488 

-480 

.2721803 

.029692 

.21034 

520 

2.879916 

-460 

.3464673 

3.582633E-02 

.2410144 

501. 

5974 

2.943744 

-440 

.4415205 

4 .529992E-02 

.2785504 

483. 

3218 

3.009938 

-420 

.5573696 

6.047795E-02 

.3231911 

465. 

1881 

3.078418 

-400 

.6872897 

.0856887 

.3764825 

447. 

2136 

3 . 149025 

-380 

.8146549 

.1287804 

.4404623 

429. 

4182 

3.221497 

-360 

.9159235 

.2030898 

.5177304 

411. 

8252 

3.29544 

-340 

.9745625 

.  3272475 

.6108596 

394. 

4617 

3. 370274 

-320 

.9958122 

.5141411 

.7233319 

377. 

3593 

3.445185 

-300 

.9997211 

.7388175 

.8597726 

360. 

5551 

3.519032 

-280 

.9999951 

.9177597 

1.025158 

344. 

093 

3. 590261 

-260 

1 

.9896971 

1.224529 

328. 

0244 

3.656775 

-240 

1 

.999714 

1.463644 

312. 

41 

3.71578 

-220 

1 

.9999992 

1.74981 

297. 

3214 

3.763606 

-200 

1 

1 

2.085394 

282. 

8427 

3.795517 

-180 

1 

1 

2.470496 

269. 

0725 

3.805514 

-160 

1 

1 

2.901352 

256. 

.125 

3 .786234 

-140 

1 

1 

3.359555 

244  . 

,1311 

3.72898 

-120 

1 

1 

3.819509 

233  . 

,2381 

3.624058 

-100 

1 

1 

4.239317 

223, 

6068 

3.461596 

-80 

1 

1 

4.532758 

215. 

.4066 

3 . 232958 

-60 

1 

1 

4.618943 

208. 

,8061 

2.932828 

-40 

1 

i 

4.44342 

203. 

,9608 

2.5617 

-20 

1 

1 

3.97739 

200. 

,9975 

2 . 128285 

0 

1 

1 

3.332139 

200 

1.651697 

20 

1 

.9999993 

2 . 540665 

200. 

.9975 

1.  165581 

40 

.9999931 

.908713 

1.76834 

203 

.9608 

.7394171 

60 

.9933111 

.467368 

1. 110188 

208 

.8061 

. 5637183 

80 

.9998733 

.7874625 

.9039734 

215 

.4066 

.7782836 

100 

.9999999 

.9856169 

1.086971 

223 

.6068 

1.134  662 

120 

1 

.9992399 

1.254863 

233 

.2381 

1.490748 

140 

1 

.9998478 

1.308583 

244 

.1311 

1.811734 

160 

1 

.9998431 

1.29482 

256 

.125 

2.089946 

180 

1 

.9994995 

1.230946 

269 

.0725 

2.326072 

200 

1 

.9974561 

1.131703 

282 

.8427 

2.523877 

220 

1 

.9872744 

1.017684 

297 

.3214 

2.688154 

240 

.999999 

.9510013 

.9018685 

312 

.41 

2.823785 

260 

.9999747 

.8646194 

.7917129 

328 

.0244 

2.935307 

280 

.9996544 

.7242736 

.6909791 

344 

.093 

3.026728 

300 

.9973592 

.5582789 

.6011183 

360 

.5551 

3.10149 

320 

.9876337 

.4037092 

.5222269 

377 

.3593 

3.162495 

340 

.9608865 

.2816961 

.4536729 

394 

.4617 

3.212155 

360 

.9086407 

.1948741 

.3944811 

411 

.8252 

3.252463 

380 

.8301703 

.1365131 

. 3443362 

429 

.4182 

3.285058 

400 

.7333605 

9.820431E-02 

.3014853 

447 

.2136 

3. 311281 

420 

.6299387 

.0731037 

.2645045 

465 

.1881 

3.332228 

440 

.5302704 

5.646636E-02 

.2325951 

483 

.3218 

3.348799 

460 

.4408869 

4 . 522933E-02 

.205043 

501 

.5974 

3.361725 

480 

.3645446 

3.747106E-02 

.1812233 

520 

3.371603 

500 

.3013791 

3.199235E-02 

.1605954 

538 

.5165 

3. 378925 

Table  3.   An  example  of  a  lateral  range  function  values  printout 
The  heading  for  the  crosscorrelation  values  is  p(cc) 
and  the  heading  for  the  square  law  values  is  p(sl). 
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Figure  3.  A  plot  of  the  square  law  lateral  range  function  values 
that  are  listed  in  Table  3.   The  horizontal  axis  is 
encounter  horizontal  range  at  CPA  in  meters.   The 
vertical  axis  is  encounter  detection  probability. 
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data. mad  magnetic  signal  values 
for  a  lateral  range  of   40   meters 

relative  CPA  distance  magnetic  signal 

s  in  meters  Hs  in  gamma 

-927.4282 

-881. 0568 

-834  .  6854 

-788.  314 

-741.9426 

-695. 5712 

-649. 1998 

-602.8283 

-556.4569 

-510.0855 

-463.7141 

-417.3427 

-370.9713 

-324.5999 

-278.2285 

-231.8571 

-185.4856 

-139.1142 

-92.74282 

-46.37141 

0 

46.37141 

92.74282 

139.  1142 

185.4856 

231.8571 

278.2285 

324.5999 

370.9713 

417. 3427 

463.7141 

510.0855 

556.4569 

602.8283 

649.  1998 

695.5712 

741.9426 

788.314 

834.6854 

881.0568 

927.4282 


-1.619177E- 

-02 

-1.834315E- 

-02 

-2 . 086144E- 

-02 

-2.381637E- 

-02 

-2.728804E- 

-02 

-.0313641 

-3.613118E- 

-02 

-4 .165412E- 

-02 

-4 .792992E- 

-02 

-5.478963E- 

-  0  2 

-6. 169327E- 

-  0  2 

-6.730653E- 

-02 

-.0686387 

-5.933568E- 

-02 

-2.652314E- 

-02 

.0540608 

.2215785 

.5163918 

.9104893 

1.184586 

1.003207 

.3837067 

-. 2252113 

-.5323231 

-.5837545 

-.5168925 

-.4202569 

-.3304064 

-.2572838 

-.2007385 

-. 1577988 

-.125285 

-.1005532 

-8. 158605E- 

-  0  2 

-6.689598E- 

-02 

-5.539922E- 

-02 

-4.630719E- 

-02 

-3 . 904362E- 

■c:1 

-3. 318428E- 

■02 

-.0284142 

-2.449729E- 

■  0  2 

Table  4.   An  example  of  a  magnetic  signal  values  printout 
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1.4 
1.3 
1.2  - 
1.1   - 

1  - 
0.9 
0.8  - 
0.7 
0.6  - 
0.5  - 
0.4 
0.3  - 
0.2  - 
0.1   - 

0 


-0  1  - 
-02  - 
-0.3  - 

-0.4 
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Figure  4.  A  plot  of  the  magnetic  signal  values  listed  in  Table  4 
The  horizontal  axis  is  relative  CPA  distance  in  meters 
The  vertical  axis  is  magnetic  signal  in  gamma. 
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program  file  name  MAD.BAS 

program  data  file  name  datal.mad 

magnetic  data  file  name  data. mag 

processing  data  file  name  data. pro 

kinematic  data  file  name  datal.kin 
combined  magnetic,  processing  &  kinematic  data  file  name  datal.mpk 
encounter  latitude  (decimal  degrees)  45 

encounter  longitude  (decimal  degrees)  -45 

encounter  variation  (decimal  degrees)  -25 
encounter  dip  angle  (decimal  degrees)  68.41608 

encounter  magnetic  field  intensity  (gamma)  59035.05 

permanent  longitudinal  moment  (oersted-cm3)  0 

permanent  transverse  moment  (oersted-cm3)  0 

permanent  vertical  moment  (oersted-cm3)  0 

target  displacement  (tons)  1800 

target  longitudinal  permeability  coefficient  13140 

target  transverse  permeability  coefficient  2880 

target  vertical  permeability  coefficient  2880 

target  longitudinal  permeability  factor  7.3 

target  transverse  permeability  factor  1.6 

target  vertical  permeability  factor  1.6 

sampling  period  (seconds)  .5 

integration  time  (seconds)  20 

adjusted  integration  time  (seconds)  20.5 

number  of  samples  per  encounter  41 

magnetometer  course  (decimal  degrees)  290 

magnetometer  speed  (knots)  180 

magnetometer  altitude  (meters)  100 

target  course  (decimal  degrees)  273 

target  speed  (knots)  10 

target  depth  (meters)  100 

magnetometer  relative  course  (decimal  degrees)  290.9828 

magnetometer  relative  speed  (knots)  170.462 

magnetometer-target  vertical  separation  (meters)  200 

target  induced  longitudinal  dipole  moment  (oersted-cm3)    1.339682E+08 
target  induced  transverse  dipole  moment  (oersted-cm3)      5. 522353E4-07 
target  induced  vertical  dipole  moment  (oersted-cm3)        1.58099E+08 
magnetic  dipole  moment  (oersted-cm3)  2.144584E+08 

dipole  moment  azimuth  (decimal  degrees)  295.4022 

dipole  moment  depression  angle  (decimal  degrees)  47.49356 

distance  between  samples  on  the  relative  track  (meters)  43.84662 
false  alarm  rate  (false  alarms  per  hour)  3 

false  alarm  probability  1.666667E-02 

magnetic  noise  (gamma)  .35 

maximum  lateral  range  (meters)  500 

lateral  range  step  (meters)  20 

number  of  lateral  range  function  values  51 


Table  5.   An  example  of  an  encounter  parameter  values  printout 
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datal.mad  lateral  range  function  values 


L 

p(cc) 

p(sl) 

H 

R 

k 

meters 

gamma 

meters 

-500 

9.718663E-02 

1.815607E- 

-02 

.1162113 

538, 

.5165 

2. 

,117317 

-480 

.1097076 

1.887764E- 

-02 

.1307452 

520 

2. 

,066198 

-460 

.1252767 

1.979502E- 

-02 

.1476683 

501 

.5974 

2. 

,017377 

-440 

.1449347 

2.098504E- 

-02 

.1677343 

483 

.3218 

1. 

,971743 

-420 

.1701671 

2.256657E- 

-02 

.1915806 

465 

.1881 

1. 

,930378 

-400 

.2031047 

2.473001E- 

-02 

.2196898 

447, 

.2136 

1. 

,89458 

-380 

.2467793 

2.779396E- 

-02 

.2543952 

429 

.4182 

1. 

,865873 

-360 

.3053632 

3.231674E- 

-02 

.2951223 

411, 

.8252 

1. 

,845999 

-340 

.3841007 

3.932704E- 

-02 

.3461031 

394 

.4617 

1. 

.836884 

-320 

.4881344 

5.082678E- 

-02 

.4085031 

377, 

.3593 

1. 

,840555 

-300 

.6185669 

7.091362E- 

-02 

.4830141 

360, 

.5551 

1. 

,859016 

-280 

.7643011 

.1082717 

.5793262 

344 

.093 

1. 

,89408 

-260 

.8946846 

.1810074 

.6993423 

328 

,0244 

1. 

,947158 

-240 

.9731281 

.3214172 

.8450639 

312, 

.41 

2. 

.019043 

-220 

.9973629 

.5584107 

1.022772 

297, 

.3214 

2. 

,109683 

-200 

.9999472 

.8329318 

1.262045 

282, 

.8427 

2. 

,217957 

-180 

.9999999 

.9794214 

1.555326 

269 

.0725 

2. 

.341425 

-160 

1 

.9996887 

1.906387 

256 

.125 

2. 

.476081 

-140 

1 

.9999999 

2.315586 

244, 

,1311 

2. 

,616057 

-120 

1 

1 

2.774995 

233, 

.2381 

2. 

.753399 

-100 

1 

1 

3.289594 

223, 

.6068 

2. 

.878014 

-80 

1 

1 

3.848771 

215. 

,4066 

2. 

,978009 

-60 

1 

1 

4.356607 

208. 

.8061 

3. 

.040665 

-40 

1 

1 

4.738054 

203. 

,9608 

3. 

,054164 

-20 

1 

1 

4.921432 

200. 

,9975 

3. 

.009893 

0 

1 

1 

4.862248 

200 

2. 

,904766 

20 

1 

1 

4.561756 

200. 

,9975 

2. 

,742697 

40 

1 

1 

4.06876 

203. 

.9608 

2. 

,534609 

60 

1 

1 

3.462592 

208, 

,8061 

2. 

.296925 

80 

1 

1 

2.876776 

215. 

.4066 

2. 

,049263 

100 

1 

.9999861 

2.341757 

223. 

,6068 

1. 

,812213 

120 

1 

.9908414 

1.862014 

233. 

,2381 

1. 

.60562 

140 

.9999632 

.8492265 

1.463262 

244. 

,1311 

1. 

.446881 

160 

.9973499 

.55795 

1.150695 

256. 

.125 

1. 

,348197 

180 

.9773154 

.3394268 

.9187382 

269. 

.0725 

1. 

,31296 

200 

.9311967 

.2232209 

.7327372 

282. 

.8427 

1, 

,334314 

220 

.8711751 

.1618173 

.5979505 

297, 

.3214 

1. 

.39856 

240 

.8071772 

.1253278 

.4948332 

312. 

.41 

1, 

,49089 

260 

.7414287 

.1006817 

.4163041 

328. 

,0244 

1, 

,599193 

280 

.6739562 

8.250537E- 

-02 

.3566572 

344. 

.093 

1, 

.714996 

300 

.6054811 

.0684997 

.3130742 

360, 

,5551 

1, 

.832886 

320 

.5376865 

5.753178E- 

-02 

.2806611 

377, 

.3593 

1, 

.949611 

340 

.4726234 

4.891159E- 

-02 

.2530364 

394. 

.4617 

2. 

.063321 

360 

.4121077 

4.213673E- 

-02 

.2280127 

411. 

,8252 

2, 

,173046 

380 

.3573765 

3.681085E- 

-02 

.2054994 

429, 

.4182 

2, 

,27835 

400 

.3090144 

3.261634E- 

■02 

.1853368 

447. 

,2136 

2. 

.379115 

420 

.2670518 

.0293013 

.1673323 

465. 

,1881 

2. 

.475412 

440 

.2311289 

2.666843E- 

■02 

.151282 

483. 

3218 

2. 

,567413 

460 

.2006615 

2.456539E- 

■02 

.1369845 

501. 

5974 

2. 

,655341 

480 

.1749714 

.0228749 

.12437 

520 

2. 

739442 

500 

.1533755 

2.150718E- 

■02 

.1132715 

538. 

5165 

2. 

819965 

Table  6.   An  example  of  a  lateral  range  function  values  printout 
The  heading  for  the  crosscorrelation  values  is  p(cc) 
and  the  heading  for  the  square  law  values  is  p(sl). 
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Figure  5.   A  plot  of  the  square  law  lateral  range  function  values 
that  are  listed  in  Table  6.   The  horizontal  axis  is 
encounter  horizontal  range  at  CPA  in  meters.   The 
vertical  axis  is  encounter  detection  probability. 
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datal.mad  magnetic  signal  values 
for  a  lateral  range  of   40  meters 

relative  CPA  distance  magnetic  signal 

s  in  meters  Hs  in  gamma 

-876.9324  6.869243E-03 

-833.0858  8.988515E-03 

-789.2391  1.181241E-02 

-745.3925  1.561336E-02 

-701.5459  2.078687E-02 

-657.6993  2.791575E-02 

-613.8527  3.787262E-02 

-570.006  5.198577E-02 

-526.1594  7.231193E-02 

-482.3128  .1020894 

-438.4662  .1464967 

-394.6196  .2139163 

-350.7729  .317984 

-306.9263  .480659 

-263.0797  .735856 

-219.2331  1.130141 

-175.3865  1.707647 

-131.5399  2.448238 

-87.69324  3.132554 

-43.84662  3.263226 

0  2.439828 

43.84662  1.023829 

87.69324  -.1487533 

131.5399  -.7045798 

175.3865  -.8055338 

219.2331  -.7078842 

263.0797  -.5631492 

306.9263  -.4307948 

350.7729  -.3256765 

394.6196  -.2465582 

438.4662  -.1881302 

482.3128  -.1451051 

526.1594  -.1132614 

570.006  -8.948089E-02 

613.8527  -7.152954E-02 

657.6993  -5 . 782279E-02 

701.5459  -4.723642E-02 

745.3925  -3 . 896839E-02 

789.2391  -3.244173E-02 

833.0858  -2.723733E-02 

876.9324  -2 . 304761E-02 


Table  7.   An  example  of  a  magnetic  signal  values  printout 
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3.5 


3  - 


2.5  - 


2  - 
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Figure  6.   A  plot  of  the  magnetic  signal  values  listed  in  Table 
7.   The  horizontal  axis  is  relative  CPA  distance  in 
meters.   The  vertical  axis  is  magnetic  signal  in 
gamma. 
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Appendix  10.  A  Program  Listing 

10  CIS  :  CLEAR 

20  PRINT  "Magnetic  Anomaly  Detection  (MAD)  Lateral  Range  Function  Program" 

30  DIMX0(70),  R0(70),  H(70)  ,  HS(70,  150),  K(70)  ,  PDCC(70)  ,  PDSL(70) 

40  PI  =  3. 141592654 #:  CON  =  2  *  PI  /  360:  KDN  =  1852  /  3600:  N$  =  "MAD.EAS" 

50  Q0  =  .2316419:  Ql  =  .31938153#:  Q2  =  -.356563782#:  Q3  =  1.781477937# 

60  Q4  =  -1.821255978#:  Q5  =  1.330274429# 

70  II  =  2.515517:  12  =  .802853:  13  =  .010328:  14  =  1.432788:  15  =  .189269:  16 

=  .001308 

80  PRINT 

90  INPUT  "generate  data  or  print  a  program  data  file  (g/p) " ;  A$ 

100  IF  A$  =  "G"  OR  A$  =  "g"  THEN  GOTO  120 

110  IF  A$  =  "P"  OR  A$  =  "p"  THEN  GOTO  3070  ELSE  GOTO  90 

120  PRINT  :  A$  =  "a" 

130  INPUT  "magnetic,  processing  &  kinematic  data  entry  by  combined  file 

(y/n)";  A$ 

140  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  2120 

150  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  160  ELSE  GOTO  130 

160  INPUT  "magnetic  data  entry  by  file/keyboard  (f/k)";  A$ 

170  IF  A$  =  "F"  OR  A$  =  "f"  THEN  GOTO  780 

180  IF  A$  =  "K"  OR  A$  =  "k"  THEN  GOTO  190  ELSE  GOTO  160 

190  INPUT  "latitude  in  decimal  degrees  (N  +)";  LAT:  lATR  =  IAT  *  CON 

200  INPUT  "longitude  in  decimal  degrees  (E  +)";  LNG:  LNGR  =  -ING  *  CON 

210  INPUT  "variation  in  decimal  degrees  (E  +)";  DEC:  DECR  =  DEC  *  CON 

220  A$  =  "a" 

230  INPUT  "input  dip  angle  (y/n)";  A$ 

240  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  260 

250  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  290  ELSE  GOTO  230 

260  INPUT  "dip  angle  in  decimal  degrees  (north  magnetic  hemisphere  +)";  DIP 

270  DIPR  =  DIP  *  CON 

280  GOTO  390 

290  IATPR  =  76  *  CON:    LNGPR  =  100  *  CON 

300  D  =  SLN(LNGR  -  LNGPR)    *  COS(IATR)  :   E  =  COS  (LNGR  -  LNGPR)    *  COS  (lATR)  :    F  = 

SIN (LATR) 

310  X  =  E:   Y  =  F 

320  GOSUB  4170 

330  T  =  T  -  (90  *  CON  -  IATPR) :  E  =  R  *  SIN(T) :  F  =  R  *  COS(T) 

340  X  =  E:  Y  =  D 

350  GOSUB  4170 

360  X  =  F:  Y  =  R 

370  GOSUB  4170 

380  DIPR  =  ATN(2  *  TAN(T)) :  DIP  =  DIPR  /  CON 

390  A$  =  "a" 

400  INPUT  "input  encounter  magnetic  field  intensity  (y/n)";  A$ 

410  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  430 

420  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  450  ELSE  GOTO  400 

430  INPUT  "encounter  magnetic  field  intensity  in  gamma";  HE 

440  GOTO  460 

450  HE  =  70000!  /  SQR(3  *  COS (DIPR)  *  COS (DIPR)  +  1) 

460  PIM  =  0:  PTM  =  0:  PVM  =  0 

470  A$  =  "a" 

480  INPUT  "input  target  permanent  dipole  moments  (y/n)";  A$ 
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490  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  510 

500  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  540  ELSE  GOTO  480 

510  INHJT  "permanent  longitudinal  moment  in  oersted-cm3  (stern-to-bow  +)";  PLM 

520  INHJT  "permanent  transverse  moment  in  oersted-cm3  (port-to-starboard  +)"; 

FTM 

530  INPUT  "permanent  vertical  moment  in  oersted-cm3  (downward  +) " ;  PVM 

540  INHJT  "target  displacement  in  tons";  WT 

550  A$  =  "a" 

560  INHJT  "input  target  permeability  coefficients  or  factors  (c/f)";  A$ 

570  IF  A$  =  "C"  OR  A$  =  "c"  THEN  GOTO  590 

580  IF  A$  =  "F"  OR  A$  =  "f"  THEN  GOTO  630  ELSE  GOTO  560 

590  INHJT  "longitudinal  permeability  coefficient  in  ogs  units";  KL 

600  INHJT  "transverse  permeability  coefficient  in  ogs  units";  KT 

610  INHJT  "vertical  permeability  coefficient  in  ogs  units";  KV 

620  FL  =  KL  /  WT:  FT  =  KT  /  WT:  FV  =  KV  /  WT:  GOTO  710 

630  INHJT  "longitudinal  displacement  factor  in  ogs  units";  FL 

640  INHJT  "transverse  displacement  factor  in  ogs  units";  FT 

650  INHJT  "vertical  displacement  factor  in  ogs  units";  FV 

660  KL  =  FL  *  WT:  Kr  =  FT*WT:  KV=FV*WT 

670  A$  =  "a" 

680  INHJT  "generate  a  magnetic  data  file  (y/n) " ;  A$ 

690  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  710 

700  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  850  ELSE  GOTO  680 

710  INHJT  "magnetic  data  file  name";  M$ 

720  ON  ERROR  GOTO  730:  GOTO  740 

730  RESUME  710 

740  OPEN  "O",  #1,  M$ 

750  WRITE  #1,  IAT,  LNG,  DEC,  DIP,  HE,  PLM,  PTM,  PVM,  WT,  KL,  KT,  KV,  FL,  FT, 

FV 

760  CLOSE 

770  GOTO  850 

780  INPUT  "magnetic  data  file  name";  M$ 

790  ON  ERROR  GOTO  800:  GOTO  810 

800  RESUME  780 

810  OPEN  "I",  #1,  M$ 

820  INHJT  #1,  IAT,  ING,  DEC,  DIP,  HE,  PLM,  PTM,  PVM,  WT,  KL,  KT,  KV,  FL,  FT, 

FV 

830  CLOSE 

840  DECR  =  DEC  *  CON:  DIPR  =  DIP  *  CON 

850  A$  =  "a" 

860  INPUT  "processing  data  entry  by  file/keyboard  (f/k)"»  A$ 

870  IF  A$  =  "F"  OR  A$  =  "f"  THEN  GOTO  1450 

880  IF  A$  =  "K"  OR  A$  =  "k"  THEN  GOTO  890  ELSE  GOTO  860 

890  A$  =  "a" 

900  INPUT  "input  sampling  period  (y/n)";  A$ 

910  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  930 

920  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  960  ELSE  GOTO  900 

930  INPUT  "sampling  period  in  seconds";  DT 

940  IF  DT  <=  0  THEN  PRINT  :  PRINT  "must  be  greater  than  zero":  PRINT  :  GOTO 

930 

950  GOTO  1090 

960  A$  =  "a" 
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970  INRJT  "input  maximum  magnetic  signal  frequency  (y/n)";  A$ 

980  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  1000 

990  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  1030  ELSE  GOTO  970 

1000  INRJT  "maximum  magnetic  signal  frequency  in  Hertz";  MAXF 

1010  IF  MAXF  <=  0  THEN  PRINT  :  PRINT  "must  be  greater  than  zero":  PRINT  :  GOTO 

1000 

1020  GOTO  1080 

1030  INPC7T  'hidnimum  target  slant  range  at  CPA  in  meters";  MINRO 

1040  IF  MINRO  <=  0  THEN  PRINT  :  PRINT  "must  be  greater  than  zero":  PRINT  : 

GOTO  1030 

1050  INRJT  "maximum  magnetometer  relative  speed  in  knots" ;  MAXVMK 

1060  MAXVM  =  MAXVMK  *  KON:  MAXF  =  2  *  MAXVM  /  MINRO 

1070  IF  MAXVM  <=  0  THEN  PRINT  :  PRINT  "must  be  greater  than  zero":  PRINT  : 

GOTO  1050 

1080  DT  =  1  /  (2  *  MAXF) :  REM  low  pass  filter,  Nyquist  sampling  rate 

1090  A$  =  "a" 

1100  INPUT  "input  integration  time  (y/n)";  A$ 

1110  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  1130 

1120  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  1210  ELSE  GOTO  1100 

1130  INPUT  "integration  time  in  seconds";  IT 

1140  IF  IT  >=  DT  THEN  GOTO  1170 

1150  PRINT  :  PRINT  "IT  =  "  +  STR$(IT)  +  "  seconds  -  minimum  =  "  +  STR$(DT)  +  " 

seconds":  PRINT 

1160  GOTO  1130 

1170  NS  =  2  *  INT(rT  /  DT  /  2)   +  1:  REM  adj  number  of  samples  per  integration 

time 

1180  IF  NS  <=  151  THEN  GOTO  1340 

1190  PRINT  :  PRINT  "IT  =  "  +  STR$(rT)  +  "  seconds  -  maximum  =  "  +  STR$(150  * 

DT)  +  "  seconds":  PRINT 

1200  GOTO  1130 

1210  INPUT  "maximum  target  slant  range  at  CPA  in  meters";  MAXR0 

1220  IF  MAXR0  <=  0  THEN  PRINT  :  PRINT  "must  be  greater  than  zero":  PRINT  : 

GOTO  1210 

1230  INPUT  •minimum  magnetometer  relative  speed  in  knots";  MINVMK 

1240  MINVM  =  MINVMK  *  KON 

1250  IF  MINVM  <=  0  THEN  PRINT  :  PRINT  "must  be  greater  than  zero":  PRINT  : 

GOTO  1230 

1260  IT  =  2  *  MAXR0  /  MINVM 

1270  IF  IT  >=  DT  THEN  GOTO  1300 

1280  PRINT  :  PRINT  "IT  =  "  +  STR$(IT)  +  "  seconds  -  iriinimum  =  "  +  STR$(DT)  +  " 

seconds":  PRINT 

1290  GOTO  1100 

1300  NS  =  2  *  INT(IT  /  DT  /  2)  +  1:  REM  adjusted  number  of  samples  per 

integration  time 

1310  IF  NS  <=  151  THEN  GOTO  1340 

1320  PRINT  :  PRINT  "IT  =  "  +  STR$(IT)  +  "  seconds  -  minimum  =  "  +  STR$(150  * 

DT)  +  "  seconds":  PRINT 

1330  GOTO  1100 

1340  A$  =  "a" 

1350  INPUT  "generate  a  processing  data  file  (y/n) " ;  A$ 

1360  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  1380 

1370  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  1510  ELSE  GOTO  1350 
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1380  INPUT  "processing  data  file  name";  P$ 

1390  ON  ERROR  GOTO  1400:  GOTO  1410 

1400  RESUME  1380 

1410  OPEN  "O",  #1,  P$ 

1420  WRITE  #1,  DT,  IT,  NS 

1430  CIOSE 

1440  GOTO  1510 

1450  INPUT  "processing  data  file  name";  P$ 

1460  ON  ERROR  GOTO  1470:  GOTO  1480 

1470  RESUME  1450 

1480  OPEN  "I",  #1,  P$ 

1490  INPUT  #1,  DT,  IT,  NS 

1500  CLOSE 

1510  A$  =  "a" 

1520  INPUT  "kinematic  data  entry  by  file/keyboard  (f/k)";  A$ 

1530  IF  A$  =  "F"  OR  A$  =  "f"  THEN  GOTO  1790 

1540  IF  A$  =  "K"  OR  A$  =  "k"  THEN  GOTO  1550  ELSE  GOTO  1510 

1550  INPUT  "magnetometer  course  in  decimal  degrees  (0  if  at  rest)";  CM 

1560  INPUT  "magnetometer  speed  in  knots";  VMK 

1570  INPUT  "magnetometer  altitude  in  meters  (below  0  is  -)";  AM 

1580  INPUT  "target  course  in  decimal  degrees";  CT 

1590  INPUT  "target  speed  in  knots";  VTK 

1600  INPUT  "target  depth  in  meters  (above  0  is  -) " ;  AT 

1610  CMR  =  CM  *  CON:  CTR  =  CT  *  CON 

1620  WXK  =  VMK  *  SIN(CMR)  -  VTK  *  SIN(CTR) 

1630  WYK  =  VMK  *  COS  (CMR)  -  VTK  *  COS  (CTR) 

1640  Z  =  AM  +  AT:  REM  vertical  separation  (-  for  magnetometer  below  target) 

1650  X  =  WXK:  Y  =  WYK:  GOSUB  4170 

1660  CR  =  T:  C  =  CR  /  CON:  REM  magnetometer  relative  course 

1670  WK  =  R:  REM  magnetometer  relative  speed 

1680  A$  =  "a" 

1690  INPUT  "generate  a  kinematic  data  file  (y/n)";  A$ 

1700  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  1720 

1710  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  1860  ELSE  GOTO  1690 

1720  INPUT  "kinematic  data  file  name";  K$ 

1730  ON  ERROR  GOTO  1740:  GOTO  1750 

1740  RESUME  1720 

1750  OPEN  "O",  #1,  K$ 

1760  WRITE  #1,  CM,  VMK,  AM,  CT,  VTK,  AT,  Z,  C,  WK 

1770  CIOSE 

1780  GOTO  1860 

1790  INPUT  "kinematic  data  file  name";  K$ 

1800  ON  ERROR  GOTO  1810:  GOTO  1820 

1810  RESUME  1790 

1820  OPEN  "I",  #1,  K$ 

1830  INPUT  #1,  CM,  VMK,  AM,  CT,  VTK,  AT,  Z,  C,  WK 

1840  CIOSE 

1850  CMR  =  CM  *  CON:  CTR  =  CT  *  CON:  CR  =  C  *  CON 

I860  MCMR  =  (CMR  -  DECR) :  REM  magnetometer  magnetic  course  in  radians 

1870  MCTR  =  (CTR  -  DECR) :  REM  target  magnetic  course  in  radians 

1880  IIM  =  KL  *  HE  *  OQS(DIPR)  *  COS  (MCTR)  :  ITM  =  -KT  *  HE  *  OOS(DIPR)  * 

SIN (MCTR) 
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1890  IVM  =  KV  *  HE  *  SIN(DIPR) 

1900  EMX  =  (PLM  +  IIM)  *  SIN(MCTR)  +  (PTM  +  ITM)  *  OOS(MCTR) 

1910  EMY  =  (PLM  +  ILM)  *  COS(MCTR)  -  (PIM  +  ITM)  *  STN(MCTR) 

1920  EMV  =  PVM  +  IVM 

1930  X  =  EMX:  Y  =  EMY 

1940  GOSUB  4170 

1950  CMLR  =  T:  REM  dipole  azimuth  relative  to  magnetic  north 

1960  X  =  EMV:  Y  =  R 

1970  GOSUB  4170 

1980  EM  =  R:  GMR  =  T:  REM  dipole  depression  angle 

1990  GML  =  CMLR  /  CON:  CM  =  OMR  /  CON 

2000  A$  =  "a" 

2010  INPUT  "generate  a  combined  magnetic,  processing  &  kinematic  data  file 

(y/n)»;  A$ 

2020  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  2040 

2030  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  2210  ELSE  GOTO  2010 

2040  INPUT  "combined  magnetic,  processing  &  kinematic  data  file  name";  E$ 

2050  ON  ERROR  GOTO  2060:  GOTO  2070 

2060  RESUME  2040 

2070  OPEN  "O",  #1,  E$ 

2080  WRITE  #1,  LAT,  LNG,  DEC,  DIP,  HE,  PLM,  PIM,  PVM,  WT,  KL,  KT,  KV,  FL,  FT, 

FV,  DT,  IT 

2090  WRITE  #1,  NS,  CM,  VMK,  AM,  CT,  VTK,  AT,  Z,  C,  WK,  ILM,  ITM,  IVM,  EM,  OML, 

OM,  M$,  P$,  K$ 

2100  CLOSE 

2110  GOTO  2210 

2120  INPUT  "combined  magnetic,  processing  &  kinematic  data  file  name";  E$ 

2130  ON  ERROR  GOTO  2140:  GOTO  2150 

2140  RESUME  2120 

2150  OPEN  "I",  #1,  E$ 

2160  INPUT  #1,  LAT,  LNG,  DEC,  DIP,  HE,  PLM,  PIM,  PVM,  WT,  KL,  KT,  KV,  FL,  FT, 

FV,  DT,  IT 

2170  INPUT  #1,  NS,  CM,  VMK,  AM,  CT,  VTK,  AT,  Z,  C,  WK,  ILM,  ITM,  IVM,  EM,  OML, 

OM,  M$,  P$,  K$ 

2180  CLOSE 

2190  DECR  =  DEC  *  CON:  DIRR  =  DIP  *  CON 

2200  CR  =  C  *  CON:  CMLR  =  OML  *  CON:  OMR  =  OM  *  CON 

2210  INPUT  "required  false  alarm  rate  in  false  alarms  per  hour";  FAR 

2220  PF  =  FAR  *  IT  /  3600:  REM  false  alarm  probability 

2230  Y  =  PF:  IF  PF  >  .5  THEN  Y  =  1  -  Y:  REM  inverse  normal  approximation 

2240  Y  =  SQR(LOG(l  /  Y  /  Y) ) 

2250  Y  =  Y  -  (II  +  Y  *  (12  +  13  *  Y) )  /  (1  +  Y  *  (14  +  Y  *  (15  +  16  *  Y) ) ) 

2260  IF  PF  <  .5  THEN  Y  =  -Y 

2270  ZP  =  -Y 

2280  CHI  =  NS  *  (1  -  2  /  9  /  NS  +  ZP  *  SQR(2  /  9  /  NS) )   3:  REM  inverse 

chi-square  approximation 

2290  INPUT  "magnetic  noise  in  gamma";  SIG 

2300  INPUT  "maximum  lateral  range  in  meters";  LRM 

2310  INPUT  "lateral  range  step  in  meters";  ST 

2320  IF  ST  <=  LRM  THEN  GOTO  2340 

2330  PRINT  :  PRINT  "maximum  step  is  "  +  STR$(LRM)  +  "  meters":  PRINT  :  GOTO 

2310 
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2340  NC  =  2  *  INT(LRM  /  ST)  +  1:  REM  number  of  lateral  range  function  values 

2350  IF  NC  <=  71  THEN  GOTO  2370 

2360  PRINT  :  PRINT  "minimum  step  is  "  +  STR$(LRM  /  35)  +  "  meters":  PRINT  : 

GOTO  2310 

2370  AIT  =  DT  *  NS:  REM  adjusted  integration  time 

2380  W  =  WK*  KON:  REM  magnetometer  relative  speed  in  meters/second 

2390  DS  =  W  *  DT:  REM  distance  between  samples  on  the  relative  track 

2400  X0  =  -(NC  -  1)  /  2  *  ST 

2410  FOR  I  =  0  TO  NC  -  1 

2420  X0(I)  =  X0 

2430  X  =  X0:  Y  =  Z 

2440  GOSUB  4170 

2450  R0  =  R:  R0(I)  =  R:  REM  target  slant  range  at  CPA  in  meters 

2460  DELR  =  T:  REM  target  depression  angle  complement  at  CPA  in  radians 

2470  IF  R0  =  0  THEN  GOTO  2800:  REM  zero  lateral  range  and  vertical  separation 

2480  EMF  =  EM  /  10  /  R0  "  3:  REM  dipole  moment  factor 

2490  CMR  =  CR  -  DECR:  REM  target  relative  magnetic  course 

2500  L  =  COS (CMR)  *  COS (CMR  -  CMLR) 

2510  M  =  COS  (DELR)  *  COS  (CMR)  *  SIN  (CMR  -  CMLR)  -  SIN  (DELR)  *  SIN  (CMR) 

2520  N  =  -SIN  (DELR)  *  COS  (CMR)  *  SIN  (CMR  -  CMLR)  -  COS  (DELR)  *  SIN(OMR) 

2530  LI  =  COS(DIPR)  *  COS  (CMR) 

2540  Ml  =  COS (DELR)  *  COS(DIPR)  *  SIN (CMR)  -  SIN (DELR)  *  SIN(DIPR) 

2550  Nl  =  -SIN (DELR)  *  OOS(DIRR)  *  SIN (CMR)  -  COS (DELR)  *  SIN(DIPR) 

2560  A2  =  2*L*L1-M*M1-N*N1:  REM  Anderson  Function  Coefficient 

2570  Al  =  3  *  (N  *  LI  +  L  *  Nl) :  REM  Anderson  Function  Coefficient 

2580  A0  =  2  *  N  *  Nl  -  L  *  LI  -  M  *  Ml:  REM  Anderson  Function  Coefficient 

2590  SUM  =  0:  HMAX  =  0:  HMIN  =  0 

2600  FOR  J  =  0  TO  NS  -  1 

2610  S  =  (J  -  (NS  -  1)  /  2)  *  DS:  BA  =  S  /  R0:  REM  Anderson  Function  Argument 

2620  AF=1/(1  +  BA*BA)   2.5:  REM  Anderson  Function  Factor 

2630  HSF  =  (A2  *  BA  *  BA  +  Al  *  BA  +  A0)  *  AF:  REM  magnetic  signal  factor 

2640  HS(I,  J)  =  DMF  *  HSF:  REM  magnetic  signal  value 

2650  IF  HS(I,  J)  >  HMAX  THEN  HMAX  =  HS(I,  J) 

2660  IF  HS(I,  J)  <  HMIN  THEN  HMIN  =  HS(I,  J) 

2670  SUM  =  SUM  +  HSF  *  HSF 

2680  NEXT  J 

2690  H(I)  =  HMAX  -  HMIN 

2700  K(I)  =  SQR(SUM) 

2710  W  =  -ZP  +  DMF  *  SQR(SUM)  /  SIG 

2720  LAM  =  EMF  *  EMF  *  SUM  /  (SIG  *  SIG)  :  AA  =  NS  +  LAM:  BB  =  1  +  LAM  /  (NS  + 

LAM) 

2730  ZN  =  -SOJR(2  *  CHI  /  BB)  +  SQR(2  *  AA  /  BB  -  1)  :  XI  =  W 

2740  GOSUB  4230 

2750  IF  Yl  >  1  THEN  Yl  =  1 

2760  PDCC(I)  =  Yl:  XI  =  ZN 

2770  GOSUB  4230 

2780  IF  Yl  >  1  THEN  Yl  =  1 

2790  PDSL(I)  =  Yl 

2800  X0  =  X0  +  ST 

2810  NEXT  I 

2820  C  =  C  /  360:  C  =  (C  -  INT(C))  *  360 

2830  IF  C  <  0  THEN  C  =  360  +  C 
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2840  OML  =  (OML  +  DEC)  /  360:  OML  =  (OML  -  INT(CML))  *  360 

2850  IF  OML  <  0  THEN  OML  =  360  +  OML 

2860  PRINT  :  A$  =  "a" 

2870  INPUT  "generate  a  program  data  file  (y/n)";  A$ 

2880  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  2900 

2890  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  3230  ELSE  GOTO  2870 

2900  INPUT  "program  data  file  name";  D$ 

2910  ON  ERROR  GOTO  2920:  GOTO  2930 

2920  RESUME  2900 

2930  OPEN  "O",  #1,  D$ 

2940  WRITE  #1,  LAT,  ING,  DEC,  DIP,  HE,  PIM,  PTM,  PVM,  WT,  KL,  KT,  KV,  FL,  FT, 

FV,  DT,  IT,  AIT 

2950  WRITE  #1,  NS,  CM,  VMK,  AM,  CT,  VTK,  AT,  IIM,  ITM,  IVM,  Z,  C,  WK,  DM,  OML, 

OM,  FAR,  PF,  SIG,  ST 

2960  WRITE  #1,  LRM,  DS,  NC,  M$,  P$,  K$,  E$ 

2970  FOR  I  =  0  TO  NC  -  1 

2980  WRITE  #1,  X0(I),  PDCC(I) ,  PDSL(I) ,  K(I) ,  H(I) ,  R0(I) 

2990  NEXT  I 

3000  FOR  I  =  0  TO  NC  -  1 

3010  FOR  J  =  0  TO  NS  -  1 

3020  WRITE  #1,  HS(I,  J) 

3030  NEXT  J 

3040  NEXT  I 

3050  CLOSE 

3060  GOTO  3230 

3070  INPUT  "program  data  file  name";  D$ 

3080  ON  ERROR  GOTO  3090:  GOTO  3100 

3090  RESUME  3070 

3100  OPEN  "I",  #1,  D$ 

3110  INPUT  #1,  LAT,  ING,  DEC,  DIP,  HE,  PIM,  PTM,  PVM,  WT,  KL,  KT,  KV,  FL,  FT, 

FV,  DT,  IT,  ATT 

3120  INPUT  #1,  NS,  CM,  VMK,  AM,  CT,  VTK,  AT,  IIM,  TTM,  IVM,  Z,  C,  WK,  DM,  OML, 

OM,  FAR,  PF,  SIG,  ST 

3130  INPUT  #1,  LRM,  DS,  NC,  M$,  P$,  K$,  E$ 

3140  FOR  I  =  0  TO  NC  -  1 

3150  INPUT  #1,  X0(I),  PDCC(I),  PDSL(I) ,  K(I) ,  H(I) ,  R0(I) 

3160  NEXT  I 

3170  FOR  I  =  0  TO  NC  -  1 

3180  FOR  J  =  0  TO  NS  -  1 

3190  INPUT  #1,  HS(I,  J) 

3200  NEXT  J 

3210  NEXT  I 

3220  CLOSE 

3230  PRINT  :  A$  =  "a" 

3240  INPUT  "print  encounter  parameter  values  (y/n) " ;  A$ 

3250  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  3270 

3260  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  3780  ELSE  GOTO  3240 

3270  LPRINT 

3280  LPRINT  "program  file  name  "  +  N$ 

3290  LPRINT  "program  data  file  name  "  +  D$ 

3300  LPRINT  "magnetic  data  file  name  "  +  M$ 

3310  LPRINT  "processing  data  file  name  "  +  P$ 
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3320 

3330 

3340 

LAT 

3350 

LNG 

3360 

DEC 

3370 

DIP 

3380 

HE 

3390 

PLM 

3400 

PIM 

3410 

PVM 

3420 

WT 

3430 

KL 

3440 

KT 

3450 

KV 

3460 

FL 

3470 

FT 

3480 

FV 

3490 

DT 

3500 

IT 

3510 

ATT 

3520 

NS 

3530 

CM 

3540 

VMK 

3550 

AM 

3560 

CT 

3570 

VTK 

3580 

AT 


LPRINT 
LPRINT 
LPRTNT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 

LPRINT 


kinematic  data  file  name  "  +  K$ 
combined  magnetic,  processing  &  kinematic  data  file  name  "  +  E$ 

encounter  latitude  (decimal  degrees)  ";  SPC(2] 

encounter  longitude  (decimal  degrees)  ";  SPC(2; 

encounter  variation  (decimal  degrees)  ";  SPC(2; 

encounter  dip  angle  (decimal  degrees)  ";  SPC(2] 

encounter  magnetic  field  intensity  (gamma)  ";  SPC(2] 

permanent  longitudinal  moment  (oersted-cm3)  ";  SPC(2; 

permanent  transverse  moment  (oersted-cm3)  ";  SPC(2] 

permanent  vertical  moment  (oersted-cm3)  ";  SPC(2] 

target  displacement  (tons)  ";  SPC(2] 

target  longitudinal  permeability  coefficient  ";  SPC(2] 

target  transverse  permeability  coefficient  ";  SPC(2] 

target  vertical  permeability  coefficient  ";  SPC(2] 

target  longitudinal  permeability  factor  ";  SPC(2] 

target  transverse  permeability  factor  ";  SPC(2; 

target  vertical  permeability  factor  " ;  SPC(2; 

sampling  period  (seconds)  ";  SPC(2] 

integration  time  (seconds)  ";  SPC(2; 

adjusted  integration  time  (seconds)  ";  SPC(2] 

number  of  samples  per  encounter  ";  SPC(2] 

•magnetometer  course  (decimal  degrees)  ";  SPC(2] 

•magnetometer  speed  (knots)  ";  SPC(2] 

•magnetometer  altitude  (meters)  ";  SPC(2] 

'target  course  (decimal  degrees)  ";  SPC(2] 

'target  speed  (knots)  ";  SPC(2; 

•target  depth  (meters)  ";  SPC(2] 
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"magnetometer  relative  course  (decimal  degrees) 

"magnetometer  relative  speed  (knots)  " 

"magnetometer-target  vertical  separation  (meters)  " 
"target  induced  longitudinal  dipole  moment  (oersted-cm3)" 

"target  induced  transverse  dipole  moment  (oersted-cm3)  " 

"target  induced  vertical  dipole  moment  (oersted-cm3)  " 

"magnetic  dipole  moment  (oersted-cm3)  " 

"dipole  moment  azimuth  (decimal  degrees)  " 

"dipole  moment  depression  angle  (decimal  degrees)  " 
"distance  between  samples  on  the  relative  track  (meters) " 

"false  alarm  rate  (false  alarms  per  hour)  " 

"false  alarm  probability  " 

"magnetic  noise  (gamma)  " 

"maximum  lateral  range  (meters)  " 

"lateral  range  step  (meters)  " 

"number  of  lateral  range  function  values  " 


•  SPC(2 
SPC(2 


3590  LPRINT 

C 

3600  LPRINT 

WK 

3610  LPRINT 

Z 

3620  LPRINT 

IIM 

3630  LPRINT 

ITM 

3640  LPRINT 

IVM 

3650  LPRINT 

EM 

3660  LPRINT 

OML 

3670  LPRINT 

CM 

3680  LPRINT 

DS 

3690  LPRINT 

FAR 

3700  LPRINT 

PF 

3710  LPRINT 

SIG 

3720  LPRINT 

LRM 

3730  LPRINT 

ST 

3740  LPRINT 

NC 

3750  FOR  I  =  0  TO  15 

3760  LPRINT 

3770  NEXT  I 

3780  PRINT  :  A$  =  "a" 

3790  INPOT  "print  lateral  range  function  values  (y/n) " ;  A$ 

3800  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  3820 

3810  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  3900  ELSE  GOTO  3790 

3820  LPRINT  D$;  "  lateral  range  function  values" 

3830  LPRINT  :  LPRINT 

3840  LPRINT  "L         p(cc)        p(sl)         H  R 

k" 
3850  LPRINT  "meters  gamma        meters" 

3860  LPRINT 

3870  FOR  I  =  0  TO  NC  -  1 

3880  LPRINT  X0(I);  TAB(10) ;  PDCC(I) ;  TAB(24) ;  PDSL(I) ;  TAB(38) ;  H(I) ;  TAB(52) ; 
R0(I);  TAB(70)  ;  K(I) 
3890  NEXT  I 
3900  PRINT  :  A$  =  "a" 

3910  INPOT  "print  magnetic  signal  values  (y/n)";  A$ 
3920  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  3940 


SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
SPC(2 
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3930  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  4120  ELSE  GOTO  3910 

3940  HUNT  :  PRINT  "identify  signal  by  encounter  lateral  range  index" 

3950  PRINT  "lateral  range  equals  the  index  times  "  +  STR$(ST)  +  "  meters" 

3960  PRINT  "index  values:  -"  +  STR$((NC  -  1)  /  2)  +  "  to  "  +  STR$((NC  -  1)  / 

2) 

3970  INPUT  "lateral  range  index";  K:  I  =  K  +  (NC  -  1)  /  2 

3980  LPRINT  :  LPRINT  :  LPRINT 

3990  LPRINT  D$;  "  magnetic  signal  values" 

4000  LPRINT  "for  a  lateral  range  of  ";  K  *  ST;  "  meters" 

4010  LPRINT 

4020  LPRINT  "relative  CPA  distance";  TAB (35) ;  "magnetic  signal" 

4030  LPRINT  "s  in  meters";  TAB(35) ;  "Hs  in  gamma" 

4040  LPRINT 

4050  FOR  J  =  0  TO  NS  -  1 

4060  LPRINT  (J  -  (NS  -  1)  /  2)  *  DS;  TAB(35) ;  HS(I,  J) 

4070  NEXT  J 

4080  PRINT  :  A$  =  "a" 

4090  INPUT  "print  magnetic  signal  values  for  a  different  lateral  range  (y/n)"; 

A$ 

4100  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  3940 

4110  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  4120  ELSE  GOTO  4090 

4120  PRINT  :  A$  =  "a" 

4130  INPUT  "continue  to  use  the  program  (y/n) " ;  A$ 

4140  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  GOTO  10 

4150  IF  A$  =  "N"  OR  A$  =  "n"  THEN  GOTO  4160  ELSE  GOTO  4130 

4160  END 

4170  R  =  SQR(X  *  X  +  Y  *  Y) :  REM  rectangular  to  polar  conversion 

4180  IF  R  =  0  THEN  T  =  0:  RETURN 

4190  IF  ABS(X  /  R)  =  1  THEN  Q  =  SGN(X)  *  (PI  /  2)  ELSE  Q  =  ATN(X  /  R  /  SQR(1  - 

X  *  X  /  R  /  R) ) 

4200  IF  ABS(Y  /  R)  =  1  THEN  T  =  (PI  /  2)  *  (1  -  SGN(Y))  ELSE  T  =  (PI  /  2)  - 

ATN(Y  /  R  /  SQR(1  -Y*Y/R/R)) 

4210  IF  Q  <  0  THEN  T  =  2  *  PI  -  T 

4220  RETURN 

4230  Yl  =  XI:  IF  XI  <  0  THEN  Yl  =  -Yl:  REM  normal  approximation 

4240  G  =  1  /  (1  +  Q0  *  Yl) 

4250  Yl  =  EXP(-Y1  *  Yl  /  2)  /  SQR(2  *  PI)  *  G  *  (Ql  +  G  *  (Q2  +  G  *  (Q3  +  G  * 

(04  +  G  *  05) ) ) ) 

4260  IF  XI  >=  0  THEN  Yl  =  1  -  Yl 

4270  RETURN 
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